Implementación de un algoritmo de solución de las ecuaciones de movimiento en dinámica molecular clásica en Unidades de Procesado Gráfico (GPU)

1 Agradecimientos Quiero aprovechar la ocasión para agradecer en especial el apoyo y constancia de mi tutor, Marco Antonio Gigosos, tanto de prácticas de empresa como del trabajo de fin de grado, en concreto la ayuda que me ha proporcionado para completar ambos proyectos, como la facilidad a acceder a ambos. Además, no podría haber llegado a realizar el trabajo de fin de grado, sin todo lo que me han enseñado y aportado otros profesores de la misma carrera en Física. Por otro lado, radicalmente opuesto, me gustaría destacar el apoyo en los campos tanto académicos como externos, pero sin los cuáles no podría haber realizado mi Grado en Física, a mis compañeros y a mi familia. 2 Resumen Este trabajo de fin de grado trata de utilizar las ecuaciones electromagnéticas correspondientes para simular la situación física correspondiente a un plasma de partículas cargadas que interaccionan entre sí. Estos programas ya existen, sin embargo, los existentes hasta ahora cuentan con una limitación procedente del hardware en la cual los bloques de la GPU no pueden sincronizarse de forma eficiente. Sin embargo, esto se soluciona en 2017 mediante la introducción en el mercado de una nueva arquitectura de GPUs. Buscaremos crear un software que aproveche las nuevas metodologías de sincronismo para que se permita hacer simulaciones con un mayor número de partículas. Se comentarán los detalles y dificultades que añade el software necesario, CUDA 9.0 o superior, con respecto a los anteriores. Palabras Clave: CUDA, Sincronismo entre Bloques, Simulación, Plasma 3 Abstract This final year dissertation will consist on using the electromagnetic equations which take part on a plasma compound of charged particles, and their way of interact between each other. Although programs like we are studying already exist, they have a hardware limitation, the GPU’s block synchronization wasn’t possible on an efective way in older devices, but since 2017 it is already solved with the introduction of new GPU’s architecture. This fact allows the simulations to be able to increase significantly the number of particles included. We will comment the details and difficulties which the new software adds concerning older ones. Keywords: CUDA, Blocks Synchronization, Simulation, Plasm 4 5 Contents 1 Introducción 8 2 Funcionamiento GPU 11 3 Equipo de Simulación 18 4 Física de la Simulación 20 4.1 Aproximaciones.......................... 20 4.2 Unidades de Simulación . . . . . . . . . . . . . . . . . . . . . 22 4.3 SituaciónInicial ......................... 24 5 Desarrollo Matemático 27 5.1 Resolución de las Ecuaciones . . . . . . . . . . . . . . . . . . 27 5.2 Situación Equilibrio . . . . . . . . . . . . . . . . . . . . . . . 29 6 Funcionamiento del Programa 30 6.1 Definición de Elementos . . . . . . . . . . . . . . . . . . . . . 30 6.2 Movimientos de memoria . . . . . . . . . . . . . . . . . . . . 34 6.3 CálculodeFuerzas........................ 39 6.4 Calculo de Posiciones . . . . . . . . . . . . . . . . . . . . . . 41 6.5 Temperatura del Sistema . . . . . . . . . . . . . . . . . . . . 42 6.6 Tiempos de Ejecución . . . . . . . . . . . . . . . . . . . . . . 44 7 Resultados 46 7.1 Primeros Pasos de Tiempo . . . . . . . . . . . . . . . . . . . 47 7.2 Temperatura del Sistema . . . . . . . . . . . . . . . . . . . . 53 7.3 Balas................................ 59 7.4 Situación de Equilibrio . . . . . . . . . . . . . . . . . . . . . 63 8 Conclusiones 67 9 Tareas Propuestas 68 6 10 Bibliografía 69 7 1 Introducción La física computacional es una técnica que ya se ha establecido como herramienta básica en muchos campos de investigación. Este tipo de experimentos que utiliza exclusivamente expresiones teóricas y una o varias maquinas que las procesen, se ha ido mejorando con los años. Poco a poco nacen herramientas informáticas que nos facilitan más el trabajo y nos permiten realizar códigos más complejos. Además, los avances en miniaturizar el hardware y hacer más eficientes los componentes, han resultado en tener mayor memoria que llenar y mayor capacidad de cálculo. Los nuevos sistemas de hardware nos permiten nuevas funciones, y nosotros las aprovechamos a la hora de realizar nuevos programas, perfilándolos para recrear mejor sistemas físicos. La simulación de partículas de plasmas tiene una gran relevancia, nos los encontramos en lugares muy diversos, como en la descarga de gases, en plasmas de astrofísica, o en situaciones concretas como en el tokamak (investigación de la fusión nuclear). Las simulaciones nos sirven para analizar características sobre las cuales no podemos extraer datos, o confirmar que las leyes teóricas se cumplen por medio de la medida experimental. Por ejemplo, podemos calcular la estadística de la velocidad las partículas en el equilibrio, el campo eléctrico medio en un punto del plasma, la cantidad de iones que se formarían en el equilibrio... Toda la utilidad que posee la simulación de partículas y el hecho de que sea necesario practicar competencias en C/C++ y aprender nuevas características de la programación en GPU (Unidad de Procesado Gráfico), que a continuación explicaremos, han sido mi motivo para realizar este TFG. Cuando hablamos de simulación, entendemos de alguna manera el resultado de los cálculos que realiza el ordenador a partir de unas ecuaciones y condiciones que introducimos. Normalmente, el ordenador ejecuta el programa en su procesador o CPU (Unidad Central de Procesamiento), de manera que realiza los cálculos de uno en uno. Hasta que no acaba uno, no empieza el siguiente, a lo que le llamamos computación secuen- 8 cial. A pesar de que las CPUs modernas se denominan multinúcleo, de manera que pueden tener hasta 2, 4, 6 o 8 hilos de cálculo activos, dependiendo de los núcleos reales o virtuales que tenga, son una cantidad de procesos simultáneos muy pequeña. En ciertas ocasiones, a la hora de realizar los cálculos, son necesarios ciertos valores anteriores propios de la simulación. Por ejemplo, a la hora de realizar los cálculos para la dinámica de una partícula necesitamos saber las posiciones de todas las partículas en el instante anterior de cálculo. Es lo que nos ocurrirá a nosotros, ya que para ejecutar la iteración n+1 necesitaremos haber finalizado la iteración n y que sus datos estén debidamente guardados. A este tipo de computación que necesita se ejecuta un paso en cuanto otro anterior acaba se denomina computación secuencial. Sin embargo, en lo que nos atañe en la simulación de partículas, en cada instante tendremos que calcular la fuerza que cada partícula ejerce sobre todas las demás partículas. Si queremos realizar esta simulación en una CPU convencional tendríamos que calcular las fuerzas de todas las partículas sobre una y después pasar a la siguiente, hasta acabar con todas ellas. Es posible realizar estos cálculos, pero es muy ineficiente, debido a que hay maneras mejores de abordar el problema. En contraposición a la computación secuencial, si tenemos una cantidad de datos conocidos y vamos a realizar distintos cálculos con ellos (¡sin cambiarlos en el proceso!) podemos ejecutar cada uno de los tratamientos de datos por separado. En el caso del cálculo de fuerzas en la simulación de partículas estamos en este caso, pues si tengo un plasma de 200 partículas, el cálculo de la fuerza que ejercen las 200 partículas sobre la número 34 es totalmente independiente a la que ejercen sobre la 122. Si les doy a dos ordenadores totalmente separados la información de la posición de todas las partículas y la situación física completa, y hacemos que cada uno calcule la fuerza sobre una partícula distinta, sería totalmente equivalente a calcular en el mismo ordenador la posición de las dos partículas, solo que invirtiendo la mitad de tiempo (despreciando el tiempo de juntar los datos al final de ambos ordenadores). En esto consiste la esencia de la 9 Para realizar el cálculo de las fuerzas vamos a necesitar realizar más de una vez sincronismos globales (de toda la malla de hilos que invocamos), ya describiremos en que puntos concretos los usaremos más adelante, en la descripción del programa. Sin considerar las nuevas versiones de CUDA tendríamos dos posibilidades: •Limitar su uso a un bloque para la simulación. Existe un mecanismo por comando de sincronismo lo suficientemente potente y rápido para permitirnos sincronizar todos los hilos de un bloque. El punto negativo es que hay una cantidad máxima de hilos que se pueden declarar en un bloque, además que la memoria compartida de acceso rápido, la memoria "Shared", es bastante pequeña, por lo que limita la cantidad de partículas que podemos ejecutar en un solo bloque. •Usar más de un bloque, pero los mecanismos de sincronización globales que existen son muy lentos, es necesario devolver el control del programa a la CPU para que lo realice y volver a ejecutar nuevo código en la CUDA. Sus ventajas son que admitiría más partículas, ya que cuentas con la memoria "Shared" de varios bloques y con los hilos que se puedan declarar en todos los bloques. Pero por otro lado se vuelve inviable de realizar por los largos tiempos de computación que produce el método de sincronismo. Pero todo cambia con la llegada de CUDA 9.0 en 2017, es un software introducido después de un cambio en hardware muy relevante. En 2016 con la entrada de la gamma de tarjetas gráficas GeForce RTX 10 hubo un cambio de arquitectura en ellas, llegando a la arquitectura Pascal. Fue un avance que permitió al software crear una malla de hilos, es decir, podemos englobar todos los hilos que invoquemos dentro de una misma malla. En la malla podemos introducir hilos de distintos bloques, incluso definir distintas mallas que se entrecrucen con los bloques. El avance que nosotros aprovecharemos es que mientras han conseguido definir una malla global de hilos, han conseguido crear una manera de sincronizar todos los hilos de la malla de manera eficiente. Utilizando un comando en código de forma análoga a como antes usabamos "syncthread()" 16 podremos realizar ejecuciones que comprendan más de un bloque y, en los puntos que nos interese, sincronizar todos nuestros procesos paralelos. Con la última versión se introducen muchas más opciones de manejo de hilos, pudiendo definir varias mallas, dividiendolas en las que nos sea necesarias a lo largo del programa, uniendolas, parando la que nos interese... Añaden mucha libertad en cuanto al manejo del hilo que procesa en sí y no solo al proceso que ejecuta. Sin embargo, a nosotros no nos hace falta todo eso, simplemente necesitamos sincronizar los hilos de toda la malla, y eso lo podemos hacer gracias a la CUDA9.0TM . 17 3 Equipo de Simulación En este apartado vamos a comentar cual será el equipo que vamos a usar, y que requerimientos o características se necesitarían para ejecutar determinados comandos obligatorios de CUDA9.0T M y superiores. A continuación, mostraremos los datos de las dos GPUs que hemos usado, las dos son idénticas, siendo ambas tarjetas gráficas GTX NVIDIA 1060 de 6GB. Para ver estos datos, el propio software de CUDA contiene unos ejemplos, entre los cuales se encuentra un programa que lee la tarjeta de procesado gráfico conectada al dispositivo. Este programa se llama "deviceQuery", resultando la tabla de la figura 2, en la que se han resaltado los datos a comentar. Figure 3: DeviceQuery - NVIDIA GeForce GTX 1060 •Recuadro 1: como vemos, se trata de la memoria "Global", la más grande de la GPU, todas las tarjetas gráficas actuales tienen como capacidad varios gigabites, en nuestro caso serán 6 Gb. Como comentamos antes, no deberíamos tener problemas de capacidad en este tipo de memoria. •Recuadro 2: aquí vemos los núcleos que tiene la gráfica. 18 •Recuadro 3: consta del tamaño de la memoria "local", la cual, aunque es pequeña, no nos dará muchos problemas de capacidad. Está relacionado con uno de los tipos de cache de la GPU. •Recuadro 4: podemos ver la asignación máxima por hardware de las memorias "Constant", "Shared" y la máxima capacidad de la memoria de registro. Suele ser la misma en todos los dispositivos, pero no está de más comprobarlo. •Recuadro 5: nos especifica el máximo número de hilos que podemos invocar, según el hardware por bloque y por multiprocesador. •Recuadro 6: cuando invocamos hilos lo podemos hacer de forma que cada uno de ellos forma parte de una celda de una malla, de modo que su índice puede depender de un dato tridimensional, este recuadro muestra las dimensiones máximas de la malla en las tres dimensiones posibles. Nosotros, por simplicidad y por no necesitar más, usaremos una malla unidimensional, todos los hilos en el "eje x". •Recuadro 7: aunque lo dejemos para el final, este es el primer sitio en el que nos tenemos que fijar para saber si podemos ejecutar el sincronismo entre bloques. Si en estas opciones nos resulta un "No", no podremos ejecutar los comandos adecuados para permitir el tipo de sincronismo que abarca más de un bloque de forma eficiente, y por lo tanto, no podremos realizar el tipo de programas descritos a continuación. No vamos a centrarnos en tanto detalle en las capacidades de la CPU porque cualquier comercial o cualquier CPU que sea compatible con una placa base que soporte una GPU moderna es capaz de ejecutar los procesos que vamos a realizar. Esto se debe a que éstos son muy sencillos, tales como imprimir datos en un archivo exterior, dar formato a variables o invocar a la GPU. Después de saber analizar nuestro dispositivo, ya podemos explicar en qué consiste la ejecución del programa en sí. 19 4 Física de la Simulación 4.1 Aproximaciones Como ya hemos mencionado previamente, vamos a aprovechar las ventajas de la computación en paralelo para simular un plasma de partículas. Todas tendrán la misma carga tanto en cantidad como en tipo de carga. Las aproximaciones que haremos serán: •La situación física consistirá en una celda cúbica de lado L, en la que tendremos las partículas de la simulación. Se supondrán condiciones periódicas en el espacio, es decir, si una partícula sale con cierta velocidad por un lado de la celda, entrará con la misma velocidad y por la cara opuesta. Es decir, supongamos que sale de la celda en la coordenada (L, L/2, L/2), la partícula volvería a entrar en la celda original por (0, L/2, L/2), justo por el mismo punto de la cara opuesta por la que sale. •Solo contaremos las interacciones de las partículas que estén como muy lejos a la mitad de lo que mide la celda. Si la partícula se encuentra en un lateral, en principio se acaba la celda, pero emularemos otra vez una celda situada en ese lugar que sea equivalente a la nuestra. Es decir, supongamos una partícula situada en las coordenadas (L- 0.1, 0, 0), una partícula situada en (0.1, 0, 0) no ejercería ninguna fuerza sobre ella, ya que está muy lejos. Sin embargo, una "partícula virtual" correspondiente a la segunda situada fuera de la celda original en las coordenadas (L + 0,1, 0, 0), sí que ejercerá fuerza sobre la primera partícula. Para contabilizar esa fuerza utilizaremos las coordenadas de la partícula que no ejerce la fuerza directamente. •La única interacción que se considerara será la de Coulomb, además, debido a la equidad de la carga de todas las partículas, esta interacción será siempre repulsiva y dependerá exclusivamente de la distancia entre ellas. ~ F=1 40π q2 r3~r (1) 20 •El tamaño de las partículas será despreciable frente al resto de distancias típicas de la simulación. Consideraremos para su dinámica el centro de masas, localizado en unas coordenadas concretas y únicas. El tamaño de las partículas, d, solo afectará a la acotación del campo eléctrico en sus proximidades inmediatas, ya que no tendremos en cuenta el choque entre partículas. A efectos prácticos, de darse el caso de choque directo entre ellas, se atravesarían como dos nubes cargadas positivamente. Sin embargo, es difícil que suceda, ya que las fuerza repulsiva tenderá a que no ocurra. La fuerza que notarán las partículas será: Si r < d, ~ F=−q2 4π0d2~r (2) Si r > d, ~ F=−q2 4π0r2~r (3) La aproximación que estamos usando hace que la fuerza deje de ser conservativa. Sin embargo, la cantidad de veces que se usará la fórmula (2) será totalmente despreciable, ya que es muy poco probable que las partículas lleguen a acercarse tanto entre sí. La distancia d será mucho menor que las distancias típicas del problema, en principio usaremos d= 10−4Unidades de Simulación de Distancia. Para asignar unas condiciones iniciales, consideraremos que todas las partículas tienen el módulo de la velocidad inicial igual, sin embargo, su dirección es fijada de forma aleatoria. Tras el paso del tiempo, las partículas irán intercambiando paulatinamente entre ellas energía cinética y potencial, obteniendo finalmente una curva de forma Maxwell Boltzmann en el histograma de su distribución de velocidades. Comprobaremos esta afirmación, además de ir viendo como llegamos a esta curva con el paso del tiempo. Podríamos haber realizado simulaciones con partículas positivas y negativas, con lo que implicaría introducir potenciales regularizados. No hemos buscado tanto la complejidad en la naturaleza del potencial, por lo que con partículas de un mismo tipo será suficiente. 21 4.2 Unidades de Simulación Para nuestra simulación, no vamos a usar unidades del sistema internacional ni similares, ya que no tienen sentido. Por un lado, no estamos resolviendo un problema totalmente real, y por el otro, la escala de unidades que nos propone el sistema internacional no nos resulta cómodo en absoluto. Así, trabajaremos con nuestro propio sistema y, de ser necesario, ya se realizará un cambio de escalas adecuado. Las magnitudes que vamos a fijar son: •Distancia: nuestra unidad de distancia típica será la distancia media entre partículas. Es decir, en cada simulación puede tener un valor distinto. •Velocidad: fijaremos la unidad de velocidad de simulación en el módulo de la velocidad inicial de nuestras partículas, ya que todas partirán con la misma velocidad. •Masa: la masa de nuestras partículas será la unidad. •Parámetro de interacción (Γ): este parámetro nos ayudará a definir la relación entre la energía potencial y la energía cinética de nuestro sistema. La energía cinética de cada partícula depende de la velocidad y la masa, siendo ambas la unidad, la energía cinética será también del orden de la unidad. Por otro lado, la energía potencial depende del parámetro de interacción y la unidad de distancia. Sabiendo que esta segunda variable es también la unidad, se puede decir que depende solo de Γ. Al estar en un plasma de partículas libres, la dinámica tiene que ser dirigida fundamentalmente por la energía cinética, siendo ella de orden superior, por lo que Γtiene que ser de orden inferior a la unidad. La expresión de Γes: Γ = q2 4π0r0 1 2kT (4) siendo r0la distancia media entre partículas. Para distintas situaciones este valor puede ser radicalmente distinto, por ejemplo, para situaciones de descarga de gases a unos cientos 22 de grados celsius, el valor de Γserá del orden de las centésimas de unidad. Si la densidad de partículas es baja, este valor tenderá a tomar valores a su vez pequeños, por ejemplo, en el tokamak disminuye hasta valores del orden de 10−6. Esto producirá un aumento en el tiempo que requerirá la simulación para alcanzar la estabilización que buscamos. En nuestro programa usaremos valores de Γtales que: 0.1>Γ>0.001 (5) Dependiendo de la intensidad de la interacción que busquemos iremos variando su valor, sin embargo manteniendo el valor en un rango correspondiente al que tendría un gas en descarga. El tiempo de la simulación, aparte del valor de la intensidad de interacción, está relacionado de forma directa con el paso de tiempo que elijamos. Si aumentamos la interacción, habrá que hacer que el paso de tiempo sea más fino de forma inversamente proporcional para mantener la calidad. Más adelante se detallará el motivo en el Desarrollo Matemático. 23 4.3 Situación Inicial Las condiciones iniciales, como hemos comentado, serán simples. Por un lado para la posición, elegiremos lugares aleatorios dentro de la celda, y por el otro, en la distribución de velocidades crearemos una delta en módulo de velocidad igual a uno, siendo la dirección de movimiento aleatoria. En los primeros instantes, debido a que las posiciones son aleatorias y las velocidades son forzadas, no tenemos porque estar en una situación de equilibrio, de hecho, es obvio que no podemos estar en ella, pues el espectro en la distribución de velocidades será distinto ya que tiene que haber interacciones entre las partículas. Si partimos de una situación inicial de no equilibrio van a suceder intercambios de energía tanto cinética como potencial. Si además la distribución de las posiciones iniciales no es la adecuada, habrá intercambio entre la energía cinética y potencial, lo que provocará un enfriamiento o un calentamiento global del sistema. En una simulación con partículas positivas y negativas, es decir, un plasma real, lo intuitivo puede ser pensar que el gas siempre se calienta, pero nada más lejos de la verdad, dependerá totalmente de la disposición inicial de las partículas. Al distribuir las partículas de distintos tipos de forma aleatoria, estas van a tender a estar repartidas de forma homogénea, es decir, tendrán energía potencial cerca de cero. En el equilibrio, la energía potencial debería ser negativa, por lo que se tendrá que producir un aumento en la energía cinética para que la potencial disminuya. 24 Sin embargo, si en la situación inicial se crea, tanto de forma casual como forzada, una distribución en la cuál la energía potencial es menor que la que tendría el sistema en equilibrio, el sistema necesitaría ganar energía potencial y el efecto sería una disminución de temperatura. De este hecho, ya se han realizado estudios previos, de los que se extraen resultados como: Figure 4: Estabilización energías desde situaciones iniciales - Classical molecular dynamics simulations of hydrogen plasmas and development of an analytical statical model for computational validity assessment. En la figura (4), procedente de artículos de este campo [3], los puntos verdes representan los estados iniciales. Las rectas rojas son diferentes estados en los que se mantiene constante la energía del sistema, por lo que si empezamos en uno de los puntos verdes solo podemos avanzar por ellos. Los puntos azules son posibles estados de equilibrio a los que llegamos en la estabilización de las partículas. Dependiendo si el sistema comienza con mayor o menor energía potencial de la que tendría en el equilibrio, el sistema se calentará o se enfriará respectivamente. 25 El primer y más importante objetivo de la GPU será calcular la fuerza que cada partícula ejerce sobre todas las demás. Para ello, cada hilo calculará la fuerza que todas las demás ejercen sobre la partícula a la que está vinculada. Por lo tanto es necesario que en algún momento del código todos los hilos sean capaces de leer las posiciones de todas las demás partículas. Podríamos hacer que lean directamente de la memoria "Global", sin embargo es tremendamente lenta. Lo que haremos será hacer que los hilos de un bloque, cada uno realizando solo un intercambio de memoria, guarden una parte de las posiciones de las partículas de la memoria "Global" en la variable auxiliar de posición situado en "Shared". De esta manera hemos conseguido que los hilos de ese bloque puedan acceder de manera rápida a las posiciones de todas las partículas que caben en la memoria auxiliar. Cuando estos hilos acaben de usar las posiciones haremos que se almacenen la información de nuevas partículas, y así hasta que recorramos todas. El proceso resultará en que todos los hilos han podido calcular la fuerza de todas otras partículas sobre la que tienen vinculada. Esta fuerza la irán almacenando en la variable de registro de fuerza para acabar volcándola en la memoria "Global". Para acceder a cada tipo de memoria declararemos distintos índices: •i absoluto ("iabs" en el código) : será el índice que nos indique el hilo que estamos usando en la malla global, de manera que recorra desde 0 hasta el número de hilos invocados menos uno. Será necesario para usarlo en la memoria "Global" ya que definiremos arrays de tantos elementos como partículas hay en la simulación. •i relativo ("i" en el código) : este índice nos muestra la posición de nuestro hilo en el bloque. Cada bloque tiene el mismo número de hilos invocados en él, de manera que el índice irá de 0 a el número de hilos por bloque menos uno. Es necesario para utilizarlo en la memoria "Shared", pues todos los hilos del bloque van a acceder a ella y cada uno tiene que hacerlo a un lugar específico. 32 •j absoluto ("jabs" en el código) : índice que va avanzando en el programa útil para identificar las partículas que copiamos a la memoria "Shared". Para cada bloque e hilo el índice j absoluto tendrá un valor distinto. •j relativo ("j" en el código) : después de copiar los datos de las partículas para interaccionar en la memoria auxiliar tendremos que calcular la fuerza que ejercen sobre las que están en la memoria de registro de cada hilo del bloque. Para el cálculo recorreremos todas las partículas almacenadas en "Shared" con el índice j. Los índices están organizados de tal manera que vamos a calcular la interacción de j sobre i. Haremos un barrido en las j en el programa y cada hilo se encarga de una i, de manera que obtenemos todas las interacciones. Si en programa invocamos la memoria el elemento 3 del array guardado en "Shared", todos los hilos en los que se invoca el comando accederan a la memoria 3 de su bloque, en otras palabras, si son de distinto bloque accederán a distintos puntos de memoria, existe una degeneración en los índices i relativo y j relativo. Esta degeneración es necesaria para acceder a la memoria "Shared", sin embargo, se puede romper utilizando los índices i absoluto y j absoluto. Es necesario romper la degeneración para acceder a la memoria "Global", ya que es común a todas. Sin embargo, al utilizar la memoria de registro o "Local" no necesitamos utilizar índices, cada hilo se autogestiona su propia memoria. Si nombro una variable de registro cada hilo accederá a la suya. La ventaja mencionada que nos proporciona 9.0 de crear mallas la vamos a necesitar, pues incluiremos todos los hilos de todos los bloques que se ejecuten en el programa en la misma malla. Ésta contendrá una cantidad de hilos como partículas tiene nuestra simulación. Figure 6: Definición de la malla en la que trabajaremos 33 6.2 Movimientos de memoria Los datos de todas las partículas se encontrarán en la memoria "Global" son muy lentas, por lo que tendremos que transferirlos a la memoria "Shared" para su uso. En la memoria de posiciones auxiliares de cada bloque vamos a guardar la información de tantas partículas como hilos se ejecutan en cada bloque. Queremos recorrer todas las partículas, por lo que, en principio, nos hará falta tantos movimientos de memoria como bloques usemos en el programa. Al principio del programa, encargaremos a cada hilo guardar los datos de la posición de su partícula vinculada a la memoria local del propio hilo. Esta variable no cambiará en ningún momento del cálculo de fuerzas. Inmediatamente después, inicializaremos a cero todas las fuerzas y copiaremos la primera ronda de partículas a la memoria auxiliar de cada bloque. La primera ronda consistirá en copiar a la memoria "Shared" de cada bloque las posiciones que están en las memorias de registro de cada bloque. En la segunda ronda copiaremos los datos de partículas que pertenecen al siguiente bloque adyacente, cuyo índice de bloque es una unidad superior al que se van a copiar las memorias. En el tercero, serán las partículas de dos unidades superiores, y así hasta que recorramos todos los bloques. En la imagen (7) se muestra el desarrollo de las interacciones entre bloques de forma completa, donde, los números simbolizan en que paso del cambio de memorias de "Global" a "Shared", es decir, primero se calcularán las interacciones de la diagonal, e iremos moviéndonos en las columnas. Sin embargo, no es necesario todos los movimientos de memoria que aparecen en la imagen, ya que aprovecharemos a calcular la fuerza del Bloque 1 sobre el Bloque 2 a la vez que calculamos la interacción del Bloque 2 sobre el Bloque 1. Reduciremos el número de iteraciones en el bucle correspondiente a interaccionar a la mitad más uno. 34 Figure 7: Máxima cantidad interacciones de los bloques Figure 8: Bloques que interaccionan de igual manera Para realizar el cálculo de ambas fuerzas de forma simultánea entrará en juego la variable auxiliar de fuerzas almacenada en "Shared" que aún no hemos mencionado. En esa variable calcularemos la fuerza que se ejerce sobre la partícula j relativa por parte de todas las partículas pertenecientes a la variable de registro de cada hilo en el bloque en el que se encuentran. Cada hilo vinculado a su partícula tendrá como función calcular la fuerza que es ejercida sobre ella misma. Además, le introduciremos la nueva tarea de calcular la fuerza que su partícula ejerce sobre cada partícula j relativa. Como los datos de la memoria j relativa están guardados en "Shared" todos los hilos podrán leer su posición y escribir sobre su fuerza. 35 En cada paso de movimientos de memoria estaremos pasando a guardar distintas partículas en la memoria "Shared", por lo que volcaremos la variable de fuerza auxiliar, distinguida en cada bloque por un j relativo, a la memoria "Global", pasando al índice j absoluto. Este proceso de cálculo de fuerzas complementario y secundario tiene su utilidad, ya que disminuye el tiempo de cálculo a la mitad. Tendremos que realizar la mitad de iteraciones en los movimientos de memoria. No añadiremos cálculos a las iteraciones que realicemos, ya que la fuerza de i relativo sobre j relativo será igual que la fuerza de j sobre i. Guardaremos en la variable de fuerzas "Local" el sumatorio de fuerzas en i, mientras que en la variable auxiliar de fuerzas en "Shared" el índice en el sumatorio será el j. Figure 9: Nuestros cálculos de interacción de blqoues Para acabar de dejarlo claro, la Figura (9) es el resumen de cómo vamos a sumar las fuerzas para calcular las totales. La línea roja representa lo que queremos conseguir, es decir, la suma de todas esas fuerzas individuales resultará en la fuerza que ejercen todas las demás partículas sobre la ligada al hilo. En cada iteración de movimientos de memoria vamos a cálcular la suma de los datos de la línea naranja de la imagen y amarilla por separado, de manera que la suma horizontal se añada a la memoria de registro, mientras que la vertical se introduzca en la memoria auxiliar. 36 El sumatorio de las líneas verticales, debido a la simetría de la tabla, como se ve gracias a las lineas azules, son equivalentes a los sumatorios de las líneas azules horizontales. Por lo que no haremos los movimientos de memoria correspondientes para calcular lo que serían las líneas horizontales azules. Simplemente, sumaremos las verticales en el lugar correspondiente. Figure 10: Movimientos iniciales descritos anteriormente y algunas inicializaciones a cero Figure 11: Movimientos iniciales descritos anteriormente Un dato importante que hay que considerar antes de mover la memoria es, si nos compensa moverla. Si la vamos a mover de una variable a otra para simplemente hacer un número pequeño de cálculos, nos cuesta el mismo tiempo que acceder a ella para copiarla a otra memoria más rápida. En los procesos anteriores sí que nos interesaba, pues el número de cálculos hechos para las fuerzas son 6402en los peores casos que hemos ejecutado, es decir, del orden de 4∗106por cada movimiento de memoria, mientras que en los casos más simples son del orden de 350 operaciones. 37 Por ejemplo, para realizar los cálculos de la velocidad a partir de la fuerza y de la posición a partir de la velocidad, que consistirá en realizar dos productos en cada hilo, lo haremos directamente en la memoria Global. En las imágenes (10) y (11) vemos como pasamos los datos de las memorias "Global" a las memorias "Shared" y "Local", siendo las de "Shared" ya dentro del bloque correspondiente. Por otro lado, en las imágenes (12) y (13) mostramos como devolvemos a memoria "Global" las variables de fuerzas, sumándolas a la número que existía previamente a esa variable. Figure 12: Devolvemos fuerzas de shared a la memoria Global en cada iteración en los bloques Figure 13: Devolvemos fuerzas de registro a la memoria Global en cada paso de tiempo 38 6.3 Cálculo de Fuerzas A pesar de todo lo que hemos comentado, es en éste apartado en el que se creará la física de la simulación. La evolución de nuestro sistema dependerá del tipo de interacción que establezcamos entre nuestros elementos, en este caso partículas. Todos los otros apartados con detalles técnicos de la simulación son necesarios y nos interesa analizarlos a fondo, sin embargo son más como un andamio si quieres pintar una fachada. La interacción entre las cargas, como ya hemos comentado, se deberá de forma exclusiva a la interacción coulombiana entre las partículas. Esta parte del código es bastante fácil de implementar una vez que tenemos los datos que las fórmulas matemáticas usan en memorias de fácil acceso. Tendremos en cuenta dos cosas: •Primero, tenemos que asegurarnos que la partícula está en la caja en la que consideramos que interacciona con cada partícula, es decir, que las partículas estén entre sí a una distancia menor de L/2. Para ello, calcularemos el vector entre la partícula del hilo y la que va a ejercer fuerza sobre ella. Si alguna de las componentes es mayor a L/2 le restaremos o sumaremos L a esa componente. Hemos conseguido arreglar las interacciones. Si una partícula i está en el borde superior de la celda y la partícula j está en la parte inferior, la distancia entre ellas es mayor que L/2, por lo que no interactuarían. Sin embargo, la partícula j situada en la parte inferior es la correspondiente a una partícula virtual en la celda inmediatamente superior a la nuestra. Esta partícula virtual sí que la tendremos que tener en cuenta para la fuerza sobre la partícula i. Lo que conseguimos al sumar L a la componente correspondiente del vector que une i con j es obtener el vector entre i y la partícula virtual correspondiente a j, por lo que podemos proseguir con el cálculo de fuerzas. 39 Figure 14: Redistribuir Partículas en el Programa •Segundo, utilizaremos una variable intermedia para almacenar la distancia entre las dos partículas y su valor al cubo. La guardaremos como variable de registro, para que sea de acceso rápido. Esta variable sirve simplemente para agilizar brevemente el cálculo. La usaremos para el cálculo de un vector secundario en el que almacenaremos de forma momentánea la fuerza que j relativo ejerce sobre i relativo. Este valor será distinto en cada hilo, pues las variables de registro como hemos dicho, aunque las llamemos igual, se ocupan en distintos sitios de memoria para cada proceso paralelo. Sumaremos a la variable de fuerza de registro, la fuerza de j sobre i. Por otro lado, añadiremos en un paso de cada hilo a la variable de fuerza almacenada en "Shared" la fuerza de i sobre j, que será el valor opuesto de j sobre i. Figure 15: Cálculo de Fuerza en el Programa 40 6.4 Calculo de Posiciones El dato de la fuerza total sobre la partícula lo almacenamos finalmente en la memoria "Global", en esta memoria calcularemos las posiciones nuevas de las partículas. Debido a que cada hilo solo tendrá que realizar un número muy límitado de operaciones, no más de 5 - 6, las realizaremos en la memoria lenta, "Global", ya que ya estaban ahí y moverlas es contraproducente. Simplemente realizaremos las siguientes operaciones: Figure 16: Cálculo de Posiciones en el Programa La fuerza que hemos calculado previamente es incompleta, ya que falta multiplicarla por el factor Γ. Lo haremos en la memoria "Global" con una simple multiplicación. Después, multiplicamos por los pasos de tiempo correspondientes y obtenemos las nuevas velocidades y posiciones de nuestras partículas. Por último, falta ver si alguna caja se ha salido de la celda. Para ello mediante la función módulo restaremos o sumaremos a cada coordenada del vector posición el tamaño de la celda para que vuelva a ella. Es análogo a olvidarnos de la partícula que ha salido y considerar una nueva que entra con la misma velocidad y por un sitio opuesto del que ha salido la primera. 41 Figure 19: Primeros pasos de la simulación Figure 20: Primeros pasos de la simulación 48 Las conclusiones que sacamos son que, a pesar de representar las mismas unidades de tiempo, no todas abandonan la delta inicial de la misma manera. El histograma correspondiente a paso de tiempo 5∗10−4se aprecia como está situado por encima del resto continuamente, es más lento a la hora de abandonar la situación inicial, y por tanto tardará más en estabilizarse. Esto se debe a que en cada paso de interacción, si el paso de tiempo es menor, las partículas recorren una distancia menor. Las fuerzas que se ejercen unas sobre otras tienden a compensar sus errores, siendo a veces más grandes y otras más pequeñas. Provocará un cambio aleatorio en la trayectoria, sin embargo, no nos preocupa ya que de primeras iba a ser aleatorio. Por otro lado, la energía del sistema aumentará de forma indefinida, pues al hacer las cuentas resulta un termino cuadrático que depende del paso de tiempo. Si el paso de tiempo es menor, este termino es menor, y disminuirá el incremento de temperatura en la simulación. Ya analizaremos los incrementos de temperatura en el siguiente apartado. Cabe destacar, aunque será explicado en un apartado posterior, que en la última imagen se puede apreciar que la velocidad cuadrática media no es la misma en todas las ejecuciones. Es decir, no todos los sistemas tienen la misma energía cinética efectiva. Esto se debe a que en nuestros programas se crean lo que llamaremos balas, no son más que partículas que adquieren una velocidad tres veces superior a la velocidad cuadrática media. Son anomalías que en algunos casos, a largo plazo, llegan a obtener cantidades nada despreciables de la energía del sistema. Provocará que el resto de partículas, no tan rápidas, tengan una velocidad cuadrática media menor. Podríamos dejar de tener en cuenta a las balas como partículas, pero es más interesante analizar por qué aparecen. El segundo análisis que vamos a realizar es si el hecho de tener una mayor cantidad de partículas provocaría que el sistema evolucione más rápido u no. La motivación para realizar esta propuesta se basa en el siguiente argumento: al aumentar el número de partículas la aproximación que hace- 49 mos será menor, pues si hay más partículas la celda será más grande, y tendremos en cuenta a más partículas para la interacción. Puesto que aumenta ese número, es intuitivo pensar que cuántas más interacciones hay, menos tardará a llegar al equilibrio. En las siguientes figuras, (21), (22) y (23), vemos la evolución descrita, y a diferencia de en el primer caso, las tres simulaciones se ejecutan de forma paralela, sin aparentemente ser ninguna más rápida que otra. Podemos deducir que el tiempo de llegar al equilibrio no depende del número de partículas, al menos de una manera considerable. Como vemos en la última imagen, la (23), el histograma correspondiente a 500 partículas es radicalmente distinta, su temperatura media no se acerca a la unidad. El origen del problema vuelven a ser las balas. A efectos prácticos es como la simulación correspondiente a 500 partículas y la de 3025 partículas estuvieran a distintas temperaturas. Vemos que son curvas situadas en lugares distintos y es más difícil comparar su forma en lugar de si las dos estuvieran dibujadas en el mismo lugar, como las correspondientes a 1860 y 3025 partículas. 50 Figure 21: Primeros Pasos de Simulación Figure 22: Primeros Pasos de Simulación 51 Figure 23: Primeros Pasos de Simulación 52 7.2 Temperatura del Sistema El factor de cambio de temperatura de la simulación irá cambiando. Si su valor es muy alto indica que la temperatura en la unidad de tiempo que se ha calculado se habría calentado más que si su valor es cercano a la unidad. A su vez, puede ser menor que la unidad, de manera que representaría un enfriamiento del sistema. Como hemos visto en la ecuación (14), la temperatura relativa a la temperatura correspondiente a la velocidad cuadrática media unidad, será simplemente la velocidad cuadrática media en el instante que queremos medirla. La energía cinética del sistema también tiene una relación biyectiva con la velocidad cuadrática media, por lo que se dará la misma situación que con la temperatura. Así pues, podremos ver el factor de cambio de temperatura como el factor de cambio de velocidad cuadrática media o el factor de cambio de energía cinética del sistema. Siendo todos los anteriores factores equivalentes tanto en sentido físico como matemático, por lo que utilizaremos los tres indistintamente. Como medimos este factor de cambio cada unidad de tiempo no dejamos al sistema suficiente libertad como para que sea muy elevado. Nos encontraremos con valores cercanos a la unidad, como podemos ver en la imagen (24). Tras la realización de todas las simulaciones, vamos a analizar de que factores depende la variación de temperatura en el sistema. La primera fuente de incremento de este factor será la magnitud de los pasos de tiempo. A continuación, compararemos los factores de enfriamiento para simulaciones con el mismo número de partículas y la misma constante de interacción, pero con pasos de tiempos distintos. Es de esperar que cuanto más pequeño sea el paso de tiempo, más cercano a la unidad sea el factor de enfriamiento, como vemos en la imagen (24). 53 Figure 24: Valores de renormalización de energía Claramente, en la figura (25) vemos cómo sucede lo que hemos descrito. Es relevante recordar que el definir pasos de tiempo muy pequeños, eleva el tiempo de computo de maneras considerables, antes de ejecutar el código hay que hacer un balance de que compensa más. Una simulación que tarde poco tiempo y una "calidad" peor, en la que se está obligado a introducir un mecanismo de enfriamiento o una simulación más fiel a la realidad en la que no haya un calentamiento artificial con la desventaja de ser mucho más lenta. Si sabemos que el plasma que estamos simulando se va a calentar hasta un equilibrio y queremos saber cómo y cuánto, se tendrá que realizar un estudio de cuanto se calienta por culpa de los pasos temporales elegidos en esa configuración para el plasma. Posteriormente elegir un paso de tiempo que nos ha proporcionado una variación de temperatura nula o despreciable en lo que el sistema se equilibra. 54 Figure 25: Variación renormalización energía para distintos pasos de tiempo El siguiente resultado corresponderá a realizar dos simulaciones, cada una con un factor de interacción entre las partículas distinto. Todos los demás parámetros de la simulación son iguales en ambas, resultando la figura (26). Cuando aumentamos la constante de interacción, aumenta de forma drástica la variable que utilizamos para enfriar el programa, es decir, la temperatura aumenta más rápidamente. Si necesitamos que nuestro parámetro de interacción sea grande, tendremos que ser consecuentes y utilizar intervalos de tiempo más finos en compensación. Por otro lado, si usamos factores Γmuy pequeños, nos da una ligera libertad de ampliar el paso de tiempo. 55 Figure 26: Valores de renormalización de energía para distintas Gammas Hay que tener cuidado, el par paso de tiempo y factor de interacción solo afecta a la velocidad de la partícula, sin embargo, para el cálculo de posiciones en cada instante solo utilizaremos el paso de tiempo. No podemos elegir pasos de tiempo que nos estropeen el cálculo de las posiciones aunque no nos estropeen las velocidades gracias al factor de interacción. Por último, en este apartado, será en el único en el que hablaremos de las ejecuciones con el programa que no enfría, si no que permite una evolución natural, pues lo que veremos en él será lo que aparece en la figura (27). Lo único que vamos a ver será que la temperatura aumenta de forma progresiva. Si que es verdad, que cuánto más grande sea el paso de tiempo la temperatura crece más rápidamente. Es posible encontrar pasos de tiempo donde el incremento de temperatura artificial, es decir culpa del paso de tiempo, sea despreciable a lo largo del programa, sin embargo, 56 Figure 27: Valores de la Temperatura sin mecanismo de enfriamiento alarga el proceso de simulación bastante. Si nuestro objetivo es realizar simulaciones con parámetros concretos y un estudio más concienzudo de los resultados, si que interesaría buscar y usar esos pasos de tiempo. Sin embargo, nosotros estamos observando que sucede cuando variamos los parámetros, nos interesa poder realizar simulaciones rápidas, por lo que introducimos el mecanismo de enfriado y elegimos pasos más gruesos, como los de la figura (27). De todos los resultados anteriores extraemos que no hay incremento de temperatura que no se deba al efecto de la discretización (que se va reduciendo a medida que el paso de tiempo se hace más pequeño), es decir, un error numérico por el hecho de haber introducido ecuaciones de diferencias. 57 La gráfica de la figura (33), correspondiente a 1860 partículas, es algo intermedia, mientras que la imagen (34), es causada por 3025 partículas, se ajusta bastante bien a la curva teórica. Es esta una de las razones, por las que necesitamos o buscamos aumentar el número de partículas, pues aumenta bastante la precisión en este tipo de datos. Figure 32: Tabla de equilibrio para 500 partículas En la gráfica (34) vemos cual es la dificultad que nos han añadido las balas, la temperatura efectiva del sistema cambia, los distintos histogramas que resultan en las distintas unidades de tiempo no se ajustan a las mismas curvas de equilibrio. Lo que nos impide comparar de forma clara lo que sucede en este tipo de simulaciones. 64 Figure 33: Tabla de equilibrio para 1860 partículas Figure 34: Tabla de equilibrio para 3025 partículas 65 Figure 35: Tabla comparativa para distintos tamaños de partículas con simulaciones de 1860 partículas con dt = 5e-4 66 8 Conclusiones Se ha conseguido realizar una simulación que, a pesar de las dificultades tratadas en este trabajo, es capaz de emular las interacciones en las que intervienen un conjunto de partículas cargadas todas de forma uniforme. Además, hemos aprovechado las ventajas que nos ofrecen las nuevas versiones de CUDA 9.0T M . Hemos implementado alguna de las mejoras ofrecidas por el software de manera que tenemos la capacidad de realizar la simulación en más de un bloque. Tener la capacidad de utilizar varios bloques para la simulación de un mismo conjunto de partículas tiene sus ventajas, como tener más potencia de cálculo y ser capaz de tener una mayor cantidad de partículas. Sin embargo, no todo son mejoras, puesto que son ejecuciones que tardan mucho más tiempo en realizar los cálculos, tanto por aumentar el número de partículas, como por los movimientos de memoria que son necesarios. Además, si utilizamos solo un bloque, los demás están libres, por lo que podemos ejecutar distintos programas en ellos, y así tener varias simulaciones activas simultáneamente. En nuestro caso, al usar hasta 30 bloques para una única tarea, no tenemos la misma capacidad. Por otro lado el TFG, motivado por el mismo alumno, ha sido de gran ayuda para el entendimiento de cómo funciona una simulación física real, de igual manera que de saciar la curiosidad en este campo. Ha conseguido fomentar una capacidad de manejo, aunque básico, de código de programación en paralelo y de su implementación en programas. La física usada, a pesar de ser sencilla, se le ha dado alguna vuelta que puede ser interesante, como el imaginar un proceso real marcado por una ecuación diferencial pasado a una ecuación en diferencias, viendo sus repercusiones. 67 9 Tareas Propuestas Para proseguir con la línea de trabajo que encamina el TFG en uno u otro aspecto, se proponen los siguientes puntos: •En primer lugar, sería destacable el hecho de seguir experimentos que hemos realizado en nuestro TFG, sin limitaciones de tiempo, de forma que se pueda simular con muchas más partículas, con pasos de tiempo minúsculos e ir probando con lo que pueda surgir para arreglar algunos de los problemas que hemos tenido, como la producción de balas. •Muchos de los procesos nombrados en las anteriores páginas pueden utilizarse para otros tipos de programas completamente distintos, de hecho, aunque no sean de simulaciones físicas, mientras sean programas que utilicen la tecnología que hemos ido mencionando, y busquen la eficiencia en la computación en paralelo. •Experimentar con el máximo de partículas que se pueden ejecutar en los bloques para sacar el mayor rendimiento a la simulación, tanto en la gráfica que hemos usado, como en otras, ya que no todas tienen las mismas prestaciones. •El uso de varios dispositivos interconectados para cualquier tipo de simulación. Hablamos del siguiente paso de mejora en la cantidad de hilos ejecutables mediante la ejecución del mismo código en distintas GPUs. De esta manera, que se introducen en la misma placa base, de manera que las gestiona la misma CPU. Se podrán comunicar de igual manera que ahora lo hacemos, con distintos comandos y teniendo algo más de cuidado. Para fijarnos si nuestra GPU sería capaz de ello, tendremos que ir al deviceQuery que ya hemos introducido y comprobar si admite el "MultiDevice Co-op Kernel Launch". 68 10 Bibliografía •CUDA Toolkit Documentation v10.1.168 (y versiones anteriores) https://docs.nvidia.com/cuda/index.html •c 2006-2010 NVIDIA Corporation. NVIDIA CUDA Programming Guide Version 3.0. •M. A. Gigosos, D. González-Herrero, N. Lara, R. Florido, A. Calisti, S. Ferri, B. Talin. Classical molecular dynamics simulations of hydrogen plasmas and development of an analytical statical model for computational validity assessment. Physical Review E 98, 033307 (2018). •Khoirudin y Jiang Shun-Liang. 2015. GPU APPLICATION IN CUDA MEMORY. Advanced Computing: An International Journal (ACIJ), Vol.6, No.2, March 2015 •Jiwei Liu. B.S. 2010. EFFICIENT SYNCHRONIZATION FOR GPGPU. Zhejiang University, Ph.D. in Electrical Engineering. •Isaac Gelado y Michael Garland. 2019. Throughput-Oriented GPU Memory Allocation. PPoPP ’19, February 16–20, Washington, DC, USA. •R.A.Fonseca, L.O.Silva, F.S.Tsung, V.K.Decyk, W.Lu, C.Ren, W.B.Mori, S.Deng, S.Lee, T.Katsouleas, and J.C.Adam. OSIRIS: A Three- Dimensional, Fully Relativistic Particle in Cell Code for Modeling Plasma Based. Springer-Verlag Berlin Heidelberg 2002. 69 •John M. Dawson. Particle Simulation of Plasmas. Department of PHysics, University of California, Los Ángeles, California. 2 Abril 1983. •Mary Thomas. COMP 605: Introduction to Parallel Computing Lecture: CUDA Shared Memory. Department of Computer Science, San Diego University. 25 de Abril de 2017. Nota: las siguientes páginas web, consistentes en foros con una comunidad que se dedica a solucionar problemas de otros usuarios, y a enseñar los conocimientos necesarios, si bien no se pueden llamar “referencias bibliográficas” en el más estricto sentido, han sido de ayuda para la realización de este trabajo de fin de grado, por lo que se merecen una cita. •https://stackoverflow.com/ •https://github.com/ •https://www.nvidia.com/es-es/about-nvidia/communities/ •https://physics.stackexchange.com/questions/55665/plasma-stability# 70