Proy2: Simulación de transferencia de calor

Se necesita una sencilla simulación por computadora que ayude a encontrar el momento de equilibro térmico de una lámina rectangular a la que se le inyecta calor constante por su borde. La lámina (en inglés, plate) corresponde a un rectángulo de dos dimensiones de un mismo material. Para efectos de la simulación, el rectángulo es dividido en \$R\$ filas y \$C\$ columnas ambas de igual alto y ancho \$h\$ como se ve en la Figura 1. Esto genera una matriz cuyas celdas son todas cuadradas, de ancho y alto \$h\$.

Cada celda de la matriz almacena una temperatura, la cual puede cambiar en el tiempo. Se usa la notación \$T_{i,j}^k\$ para indicar la temperatura de la celda ubicada en la fila \$i\$, columna \$j\$, en el instante o estado \$k\$. Después de transcurrido un tiempo \$\Delta t\$, la simulación pasará del instante \$k\$ al instante \$k+1\$, y la temperatura en la lámina habrá variado (estado). Como es sabido, la energía se transfiere de un área más caliente hacia una más fría. La nueva temperatura en la celda \$i,j\$ será \$T_{i,j}^{k+1}\$ como se ve en la parte derecha de la Figura 1, y puede estimarse a partir de su temperatura en el instante (o estado) anterior y la temperatura de sus celdas vecinas por la relación:

De acuerdo a la relación anterior, la temperatura de una celda en el instante o estado \$k+1\$ indicado por \$T_{i,j}^{k+1}\$, es el resultado de la temperatura que la celda tenía en el instante o estado anterior \$T_{i,j}^k\$ más la pérdida o ganancia de energía que la celda haya sufrido con sus inmediaciones durante ese período \$\Delta t\$. Para efecto de la simulación, las inmediaciones son las cuatro celdas vecinas en forma de cruz como se ve en la Figura 1. Esta transferencia de energía o calor está regida por:

La Figura 2 muestra cuatro instantes o estados de una simulación hipotética de una lámina de oro (difusividad térmica \$α=127\frac{mm^2}{s}\$). Para efectos de la simulación, la lámina fue dividida en 5 filas y 4 columnas, cuyas celdas son de \$h=1000mm\$ de lado, es decir, de un metro de ancho por un metro de alto.

En el estado o instante cero (\$k=0\$) la simulación carga la matriz de un archivo que indica las temperaturas iniciales de cada celda de la lámina. Es importante resaltar que los bordes de la lámina no cambian su temperatura en el tiempo, dado que es el punto donde los experimentadores "inyectan o retiran calor". Por esto los bordes se resaltan con color de fondo en la Figura 2. De esta figura puede verse que en la parte superior se inyecta calor a una temperatura constante de 10 unidades (Celcius, Farenheit, o Kelvin), y conforme se desciende en la lámina, se provee menos calor en los bordes.

En cada instante o estado, la simulación debe actualizar las celdas internas de la lámina de acuerdo al modelo físico presentado en la sección anterior. En el instante o estado \$k=1\$ habrán transcurrido \$k \Delta t = 1200s = 20min\$. Como puede verse en la Figura 2, las temperaturas en los bordes se mantienen constantes, pero las celdas internas han adquirido energía de los bordes, en especial las celdas en la parte superior. Sin embargo, la celda \$2,2\$ perdió energía pese a que está al lado de un borde de temperatura 6, dado que tres de sus vecinas estaban más frías que ella en el estado previo \$k=0\$.

En el estado \$k=2\$ habrán transcurrido \$k \Delta t = 2 \cdot 1200s = 40min\$. Como puede verse en la Figura 2, las celdas internas han incrementado lentamente su temperatura dado a que son de \$1m^2\$ cada una. Incluso la celda \$2,2\$ ha visto reflejado un incremento. En el estado \$k=3\$ que en la vida real ocurriría una hora después de que inicia el experimento, la temperatura interna sigue creciendo, pero aún no se ha equilibrado con los valores de los bordes.

Se desea que la simulación continúe hasta que se haya alcanzado el punto de equilibrio, lo cual ocurre cuando el calor se ha estabilizado en la lámina. Para esto se proveerá un parámetro épsilon (ε) a la simulación, que representa el mínimo cambio de temperatura significativo en la lámina. En cada estado \$k\$ se actualizan todas las celdas internas de la lámina. Si al menos una de las celdas internas tiene un cambio en su temperatura mayor a ε, indica que no se ha alcanzado aún el equilibrio y la simulación continúa con el siguiente estado \$k+1\$, de lo contrario se detiene y reporta los resultados de la simulación. Por ejemplo, si la simulación de la Figura 2 se corriera con un \$ε=2\$ unidades de temperatura, ésta terminaría en el estado \$k=2\$, dado que el cambio de temperatura más grande del estado \$k=1\$ a \$k=2\$ se da en la celda \$1,1\$, calculada como \$|4.51-2.74|=1.77\$, y es menor que el \$ε=2\$.

El modelo físico presentado en la sección anterior es muy sensible a los parámetros de entrada, y dependiendo de la combinación de valores puede producir resultados incorrectos. El modelo se acerca más a la realidad entre más celdas se usen para representar la lámina (filas y columnas) y más pequeños sean los cambios de tiempo (\$\Delta t\$). Sin embargo acercarse a la realidad impone más presión sobre los recursos de la máquina, lo que hace la simulación más lenta, por lo que se desea una versión paralelizada de la simulación, que pueda encontrar el punto de equilibrio térmico en el menor tiempo posible.

Se necesita que la solución concurrente sea invocada desde la línea de comandos con los siguientes argumentos:

Por ejemplo, la siguiente invocación indica que se quiere realizar todas las simulaciones indicadas en el archivo de trabajo job001.txt con 16 hilos de ejecución.

El archivo de trabajo (job file) es un archivo de texto que lista varias láminas y los parámetros de simulación que los experimentadores quieren investigar en cada una de ellas. El siguiente es un ejemplo de un archivo de trabajo job001.txt:

Cada línea del archivo de trabajo contiene una simulación independiente de las demás. Una simulación consta de los siguientes parámetros separados por espacios en blanco:

Su programa debe realizar todas las simulaciones indicadas en el archivo de trabajo. Una simulación involucra cargar la lámina como una matriz en memoria, y actualizarla a lo largo de varios estados hasta que se haya alcanzado el punto de equilibrio térmico. Una vez que esto ocurre se debe hacer un reporte a los investigadores, como se indica en la próxima sección.

El archivo de reporte provee estadísticas resultado de ejecutar cada simulación. Su nombre tiene la forma job###.tsv, donde ### es el mismo número del archivo de trabajo. Cada vez que se invoca el programa de simulación, se crea o sobrescribe el archivo de reporte correspondiente. El archivo de reporte es de texto y tiene un formato similar al archivo de trabajo:

Las líneas del archivo de reporte coinciden con las del archivo de trabajo. Cada línea del reporte contiene los mismos valores del archivo de trabajo, pero separados por tabuladores, y agrega dos resultados:

Para formatear el tiempo transcurrido en segundos, puede usar la siguiente función:

Los experimentadores están también interesados en conocer el estado de la lámina una vez que haya alcanzado el punto de equilibrio. Por eso la simulación debe crear –o sobrescribir si ya existe– un archivo binario con el estado de la matriz en el punto de equilibrio, y con el nombre plate###-k.bin donde k corresponde al número de estado donde se alcanzó el equilibrio. Este archivo tiene el mismo formato que cualquier otro archivo de lámina, y por lo tanto, podría ser usado como punto de partida de nuevas simulaciones.

Las láminas son matrices que se almacenan en archivos binarios que tienen una estructura sencilla. Los primeros ocho bytes son un entero sin signo \$R\$ que indica la cantidad de filas de la matriz. Los segundos ocho bytes son un entero sin signo \$C\$ que indica la cantidad de columnas de la matriz. A partir de los 16 bytes anteriores, continúan \$R\cdotC\$ números flotantes de doble precisión con la temperatura inicial de cada una de las celdas de la matriz, en el orden habitual de filas y columnas. Todos los valores se encuentran en little-endian.

Importante: Las matrices que los investigadores simulen pueden ser muy grandes debido a la necesidad de una alta granularidad para poder acercar la fidelidad de la simulación a la realidad. Su programa debe hacer un cuidadoso y eficiente manejo de memoria y de errores.

Se recomienda primero resolver el problema con una solución serial enfocada en corrección. El equipo debe escoger el paradigma dominante para su solución serial entre programación procedimental (recomendado) o programación orientada a objetos (OOP).

Los diseños deben guardarse en la subcarpeta design/ dentro de la carpeta para el proyecto02 en su repositorio de control de versiones. Si usa diagramas, expórtelos a imágenes vectoriales (SVG) o escalares (PNG) e incrústelas en un archivo design/readme.ext, donde la extensión depende del tipo de documento (Markdown, AsciiDoc o LaTeX). En este documento se escriben las explicaciones que ayuden a comprender la solución (los diagramas) a alto nivel.

Implementen la versión serial en C++ aplicando las buenas prácticas de programación y documentación. La aplicación debe validar las entradas, como números mal formados, fuera de rango, matrices incompletas o inválidas. El programa debe reportar estas anomalías con mensajes de error, en lugar de caerse o producir resultados incorrectos. Asegúrense de que la aplicación no genere diagnósticos del linter o sanitizers.

Se recomienda crear un programa adicional (herramienta) que compare dos archivos binarios de lámina. El programa recibe tres argumentos en línea de comandos: el nombre del primer archivo a comparar, el nombre del segundo archivo, y una sensitividad (un épsilon). Este tercer parámetro es necesario dados los problemas de precisión que presenta la aritmética flotante. La sensitividad es un número pequeño de tal forma que una diferencia menor a este número se considera despreciable y que los dos números son el mismo.

El programa comparador compararía las matrices en ambos archivos. Si tienen las mismas dimensiones, procedería a comparar celda por celda. Las celdas son temperaturas representadas en números flotantes de doble precisión. La comparación no debe hacerse directamente con el operador de igualdad (==) o desigualdad (!=), sino que por errores de precisión aritmética, podrían haber diferencias que dependen del orden en que se realicen las operaciones o por otros motivos, aunque realmente los modelos sean correctos. Por eso se puede comparar la diferencia entre cada pareja de celdas de la matriz contra un valor de sensitividad o épsilon, como se muestra a continuación:

Este programa permitirá comparar casos de prueba binarios contra una salida esperada. El código fuente de este programa debe estar en su propia carpeta para no interferir con el código de la simulación.

Una vez construida su solución serial, realicen una paralelización que incremente el desempeño usando la tecnología OpenMP. Dado que la paralelización es una optimización, apliquen el método sugerido para optimizar visto durante las lecciones.

Paso 1. Se debe medir el rendimiento de la versión serial con un archivo de trabajo grande, siguiendo el procedimiento indicado en la Tarea03 para realizar mediciones de rendimiento en un nodo esclavo del clúster Arenal con al menos tres corridas (no usar una computadora local). Anotar los resultados en una hoja de cálculo. Nota: si no se implementó una versión serial puede usar la versión paralelizada con un único hilo.

Paso 2. Hacer profiling con Callgrind para identificar las regiones de código que más impactan el consumo de CPU con un archivo de trabajo pequeño. Este paso debería ayudar al equipo a confirmar sus hipótesis de cuáles líneas de código son las candidatas a ser paralelizadas.

Paso 3. Diseñar e implementar la paralelización con OpenMP. Discutan la optimización pretendida y el método de mapeo que logre la mayor eficiencia. Consideren que las láminas también podrían ser muy altas y delgadas, o muy anchas y aplastadas. Implementar la paralelización con OpenMP en el código fuente. Recuerde que la cantidad de hilos debe poderse controlar por argumento de línea de comandos. En caso de omitirse se debe suponer la cantidad de CPU disponibles en el equipo donde corre el programa.

Paso 4. Pasar las pruebas de corrección. La solución debe pasar los casos de prueba usando el programa comparador, y no tener malas prácticas como espera activa, condiciones de carrera, bloqueos mutuos, inanición, serialización innecesaria, o ineficiencia.

Paso 5. Medir el rendimiento de la versión paralelizada con los mismos métodos que el paso 1. Calcular el incremento de velocidad (speedup) y la eficiencia. Cree un gráfico que incluya en el eje-x las dos soluciones (serial, concurrente), en el eje-y primario el incremento de velocidad, y en el eje-y secundario la eficiencia. Para la versión concurrente use tantos hilos como núcleos hayan disponibles en la máquina. Incruste su gráfico en una sección "Análisis de rendimiento" de su readme.ext con su correspondiente discusión corta (500 palabras máximo).