Propagador de ondas acústicas paralelo en 2D/3D basado en el método pseudo-espectral: implementación y análisis de rendimiento

Facultad de Inform´ atica de Barcelona Grado en Ingenier ´ ıa Inform´ atica Especialidad de Computaci´ on Propagador de ondas ac´usticas paralelo en 2D/3D basado en el m´etodo pseudo-espectral: implementaci´on y an´alisis de rendimiento Trabajo de Fin de Grado Adri´an Matas Ruiz Director: Octavio Castillo Reyes Tutor GEP: Jorge Enrique Esteban P´erez 18 de junio de 2021 Resumen El m´etodo pseudo-espectral (PSM) es una estrategia num´erica utilizada para la soluci´on de ecuaciones de derivadas parciales (PDE) en diversos y variados campos como la din´amica de fluidos, simulaci´on de ondas no lineales y modelizaci´on s´ısmica, entre otros. En el contexto de las aplicaciones geof´ısicas, los algoritmos basados en el PSM han demostrado mayor precisi´on y eficiencia en comparaci´on con otros m´etodos num´ericos como el de diferencias finitas (FD). Sin embargo, la implementaci´on eficiente del m´etodo PSM requiere de algoritmos ´optimos y de estrategias avanzadas de paralelizaci´on. En este proyecto se analiza y mejora el flujo de trabajo de un prototipo PSM para la simulaci´on de ondas ac´usticas en 2D/3D. En particular, se estudian estrategias computacionales que impacten positivamente en la eficiencia del m´etodo num´erico y alternativas que mejoren la robustez y flexibilidad del c´odigo. El prototipo, escrito en Python, es capaz de generar im´agenes geof´ısicas que representan mapas detallados del interior de la tierra. Estas im´agenes y conocimiento tienen aplicaci´on directa en diferentes sectores de los que destacan la monitorizaci´on de fuentes de energ´ıa geot´ermica, la modelizaci´on y caracterizaci´on de yacimientos de agua, y la propagaci´on de ondas en medios porosos, entre otras. 1 Resum El m`etode pseudo-espectral (PSM) es una estrategia num´erica utilitzada per la resoluci´o d’equacions de derivadas parcials (PDE) en diversos i variats camps com la din`amica de fluids, simulaci´o d’ones no lineals i modelitzaci´o s´ısmica, entre d’altres. En el context de les aplicaci´ons geof´ısiques, els algoritmes basats en el PSM han mostrat major precisi´o i eficiencia en comparaci´o amb altres m`etodes num´erics com el de difer`encies finites (FD). Per`o, l’implementaci´o eficient del m´etodoe PSM requereix d’algoritmes ´optims i d’estrategies avan¸cades de paral·lelitzaci´o. En aquest projecte s’analitza i es millora el fluxe de treball d’un prototip PSM per la simulaci´o d’ones ac´ustiques en 2D/3D. En particular, s’estudien estrategies computacionals que impacten positivament en l’eficiencia del m`etode num´eric i alternatives que millorin la robustesa i flexibilitat del codi. El prototip, escrit en Python, es capa¸c de generar imatges geof´ısiques que representen mapes detallats de l’interior de la terra. Aquestes imatges i coneixement tenen aplicaci´o directa en diferents sectors dels quals destaquen monitoritzaci´o de fonts d’energ´ıa geot´ermica, modelitzaci´o i caracteritzaci´o de jaciments d’aigua, i la propagaci´o d’ones en medis porosos, entre d’altres. 2 Abstract The pseudo-spectral method (PSM) is a numerical strategy to sol- ve partial differential equations (PDE) applied in diverse and various fields like fluid dynamics, simulation of non-linear waves, and seismic modeling, among others. In the context of geophysical applications, the algorithms based on the PSM have shown greater precision and efficiency compared with other numerical methods such as finite difference (FD). However, an efficient implementation of the PSM method requires optimal algorithms and advanced paralelization strategies. In this project we analize and improve the main work-flow of a PSM prototype to simulate acoustic waves in 2D/3D. In particular, computational strategies that positively impact efficiecy of the numeric method and alternatives to improve the robustness and flexibility of the code are studied. The prototype, written in Python, is able to generate geophysical images that represent detailed maps of the interior of the earth. These images and knowledge have a direct application in various sectors of which the monitoring of geothermal energy sources, modeling and characterization of water reservoirs, and waves propagation on porus medium, stand out among others. 3 ´ Indice 1. Introducci´on 11 1.1. Contextualizaci´on . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.2. Actores implicados . . . . . . . . . . . . . . . . . . . . . . . . 11 1.3. Descripci´on de los problemas . . . . . . . . . . . . . . . . . . . 12 1.4. Justificaci´on............................ 12 1.5. Objetivo.............................. 12 1.5.1. Subobjetivos . . . . . . . . . . . . . . . . . . . . . . . 13 1.6. Alcance .............................. 13 1.6.1. Requerimientos funcionales . . . . . . . . . . . . . . . . 13 1.6.2. Obst´aculos y riesgos . . . . . . . . . . . . . . . . . . . 13 1.7. T´erminos y conceptos . . . . . . . . . . . . . . . . . . . . . . . 14 1.7.1. M´etodo de fourier pseudo-espectral . . . . . . . . . . . 14 1.7.2. M´etodo de diferencias finitas . . . . . . . . . . . . . . . 14 2. Metodolog´ıa y rigor 14 2.1. Herramientas ........................... 15 3. Planificaci´on temporal 16 3.1. Descripci´on de las tareas . . . . . . . . . . . . . . . . . . . . . 16 3.1.1. Gesti´on del proyecto (GP) . . . . . . . . . . . . . . . . 16 3.1.2. Trabajo previo (TP) . . . . . . . . . . . . . . . . . . . 17 3.1.3. Desarrollo e implementaci´on de mejoras (DM) . . . . . 17 4 3.2. Recursos.............................. 18 3.3. Representaci´on gr´afica de la planificaci´on . . . . . . . . . . . . 20 3.4. Gesti´on del riesgo . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.5. Cambios respecto a la planificaci´on inicial . . . . . . . . . . . 21 3.5.1. Tareas........................... 21 3.5.2. Representaci´on gr´afica de la planificaci´on actual . . . . 23 4. Presupuesto 24 4.1. Identificaci´on de los costes . . . . . . . . . . . . . . . . . . . . 24 4.1.1. Personales......................... 24 4.1.2. Directos .......................... 25 4.1.3. Indirectos ......................... 26 4.1.4. Contingencia . . . . . . . . . . . . . . . . . . . . . . . 26 4.1.5. Imprevistos ........................ 27 4.1.6. Total............................ 27 4.2. Control de gesti´on . . . . . . . . . . . . . . . . . . . . . . . . 27 4.3. Cambios respecto al presupuesto inicial . . . . . . . . . . . . . 28 5. Informe de sostenibilidad 28 5.1. Dimensi´on ambiental . . . . . . . . . . . . . . . . . . . . . . . 29 5.2. Dimensi´on econ´omica . . . . . . . . . . . . . . . . . . . . . . . 29 5.3. Dimensi´on social . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.4. Matriz de sostenibilidad . . . . . . . . . . . . . . . . . . . . . 30 5 6. Aspectos legales 30 6.1. Propiedad intelectual . . . . . . . . . . . . . . . . . . . . . . . 30 6.2. Licencias.............................. 31 7. Propagador de ondas ac´usticas 2D/3D (PSM) 32 7.1. Algoritmo de migraci´on de tiempo inverso . . . . . . . . . . . 33 7.2. Estructura del c´odigo . . . . . . . . . . . . . . . . . . . . . . . 34 7.2.1. Dependencias . . . . . . . . . . . . . . . . . . . . . . . 34 7.2.2. Entrada y salida . . . . . . . . . . . . . . . . . . . . . 35 7.2.3. Preprocesado . . . . . . . . . . . . . . . . . . . . . . . 37 7.2.4. Kernel........................... 38 7.2.5. Postprocesado . . . . . . . . . . . . . . . . . . . . . . . 39 7.3. Simulaciones y resultados . . . . . . . . . . . . . . . . . . . . 40 7.3.1. Test de validaci´on: modelo impulso-respuesta . . . . . . 40 7.3.2. Precisi´on num´erica . . . . . . . . . . . . . . . . . . . . 40 7.3.3. Tiempos y consumo de memoria en un ordenador personal............................ 42 7.3.4. An´alisis de escalabilidad . . . . . . . . . . . . . . . . . 43 8. Implementaci´on de mejoras 49 8.1. Docker............................... 49 8.2. Flujo de trabajo del c´odigo . . . . . . . . . . . . . . . . . . . . 49 9. Conclusiones y trabajo a futuro 52 6 Referencias 54 Anexo 56 Anexo A: Test Escalabilidad . . . . . . . . . . . . . . . . . . . . . . 56 7 ´ Indice de figuras 1. Diagrama de Gantt creado con Gantter. Elaboraci´on propia. . 20 2. Diagrama de Gantt creado con Gantter. Elaboraci´on propia. . 23 3. params.yaml ............................ 36 4. Ejemplo de cuadr´ıcula espacio-tiempo para un modelo 2D (izquierda) y boceto de cuadr´ıcula con ppw=2 (derecha). Elaboraci´on propia ........................... 38 5. Visualizaci´on de la precisi´on num´erica. Elaboraci´on propia. . . 41 6. Test de tiempo de ejecuci´on y consumo de memoria del prototipo. 42 7. Observable 1 (2D). Elaboraci´on propia. ............. 44 8. Observable 2 (2D). debug snapshot time skip = 2 (izquierda) ydebug snapshot time skip = 6 (derecha). Elaboraci´on propia. 45 9. Observable 3 (2D). Elaboraci´on propia. ............. 46 10. Observable 1 (3D). Elaboraci´on propia. ............. 47 11. Observable 2 (3D). debug snapshot time skip = 2 (izquierda) ydebug snapshot time skip = 6 (derecha). Elaboraci´on propia. 48 12. Observable 3 (3D). Elaboraci´on propia. ............. 48 13. Ejemplo de slab decomposition con 4 cores. (a) descomposici´on en el eje Y; (b) descomposici´on en el eje X. Tomado de la p´agina de 2decomp&FFT[12]. ................ 50 14. Ejemplo de pencil decomposition con 4*3 cores. (a) X-pencil; (b) Y-pencil; (c) Z-pencil. Tomado de la p´agina de 2decomp&FFT[12]. ................................... 50 15. Ejemplo de comunicaci´on para la interpolaci´on con halo = 3. Elaboraci´on propia. ........................ 51 8 A partir de esta etapa se utilizar´a la metodolog´ıa ´agil Scrum, ya que ´esta nos permitir´a ser flexibles e ir introduciendo cambios cada cierto tiempo pudiendo as´ı reaccionar a imprevistos. Durante este periodo se definir´an ciclos de desarrollo cortos de 2 semanas aproximadamente. Al final de estos ciclos se har´a una reuni´on mediante Google Meet[3] donde se pondr´a al d´ıa de los cambios introducidos en el proyecto al director de ´este y se discutir´an los siguientes pasos a seguir. 2.1. Herramientas Una de las herramientas esenciales utilizada en este proyecto es Git[4], un sistema de control de versiones (VCS) que nos permite introducir cambios en el c´odigo pero manteniendo las versiones anteriores permitiendo as´ı volver hacia atr´as en todo momento. Esta herramienta es muy ´util junto con GitHub[5], una p´agina web donde se pueden alojar en repositorios p´ublicos o privados estas versiones y compartirlas con quien se desee. En nuestro caso existe un repositorio privado, con la idea de hacerlo p´ublico cuando el proyecto est´e finalizado, compartido con el director del proyecto para que as´ı pueda tener un buen seguimiento. Otra herramienta que se utiliza es docker[6] que es un software que automatiza el despliege de aplicaciones en contenedores. Esto permite utilizar la aplicaci´on en un entorno controlado para garantizar al usuario el uso correc- to de esta ya que todas las dependencias usar´an la versi´on correcta en este entorno. Adem´as, se utiliza Google Meet[3] para el servicio de videollamadas gratuito. Esta herramienta permite al director crear eventos con fecha y hora para poder encontrarse v´ıa online con el desarrollador y hacer seguimiento. Tambi´en permite mostrar la pantalla a la vez que se habla lo que resulta muy ´util para poder resolver dudas de una forma gr´afica 15 3. Planificaci´on temporal Con tal de llevar a cabo correctamente el trabajo es necesario hacer una buena planificaci´on temporal. En esta secci´on se muestra dicha planificaci´on. Pese a haber matriculado el TFG en febrero, se empez´o mucho antes con el trabajo. Concretamente el proyecto comenz´o el 13 de Octubre del 2020 y se tiene previsto como fecha l´ımite acabarlo el 28 de Mayo del 2021. En total el trabajo se realiza durante 227 d´ıas y se estima que la duraci´on total del trabajo son 540 horas debido a que la carga acad´emica del TFG equivale a 18 cr´editos (30 horas por cr´edito). 3.1. Descripci´on de las tareas A continuaci´on se ponen en detalle las tareas del proyecto resultantes de los objetivos expuestos anteriormente. Las tareas se distribuyen en diferentes grupos seg´un a la fase del trabajo que correspondan. 3.1.1. Gesti´on del proyecto (GP) Este apartado incluye todas las tareas no t´ecnicas relacionadas con la asignatura de Gesti´on de Proyectos (GEP). Alcance y contexto (GP1) Esta fase pretende revisar el estado del arte y poner el proyecto en contexto de la FIB y definir su alcance especificando los objetivos, subobjetivos y requisitos del trabajo. Su duraci´on se estima que sean 25 horas. Planificaci´on temporal (GP2) En esta tarea se propone una planificaci´on temporal inicial donde se definen las tareas a realizar durante el TFG y se distribuyen en el tiempo con tal de contar con un plan desde el principio. Se estima que esta tarea dure 15 horas. Presupuesto y sostenibilidad (GP3) La tarea pretende definir los costes del proyecto con un presupuesto as´ı como al escrito de un informe de sostenibilidad econ´omica, ambiental y social. Su duraci´on se estima que sean 15 horas. 16 Memoria (GP4) Como su nombre indica, consiste en escribir la memoria final del proyecto donde consta la documentaci´on de todas las fases del proyecto as´ı como la conclusi´on final. Se prevee que tenga una duraci´on de 80 horas. Presentaci´on (GP5) Una vez finalizada la escritura de la memoria, ´esta se presentar´a ante el tribunal del TFG. En esta fase se crear´a el material de la presentaci´on como su gui´on y sus ensayos. La duraci´on de la tarea es de 20 horas. Reuniones (GP6) Durante todo el proyecto se cuenta con reuniones con el director para ponerse al d´ıa y discutir cambios o mejoras. Se prevee que estas reuniones se hagan peri´odicamente cada dos semanas aproximadamente y que la duraci´on de cada una sea de 1 hora. En total contamos con 16 horas. 3.1.2. Trabajo previo (TP) Preparaci´on del entorno de trabajo (TP1) Para poder ejecutar el c´odigo correctamente es necesario contar con entorno en el equipo. Concretamente se necesita una distribuci´on de Linux que en nuestro caso es Ubuntu18.04 LTD [7] y adem´as es necesario instalar las dependencias de la aplicaci´on. Esta tarea se estima que sea de 25 horas. 3.1.3. Desarrollo e implementaci´on de mejoras (DM) Estudio de sistema de control de versiones (DM1) En esta tarea se realizar´a un estudio del funcionamiento de Git[4] yGitHub[5] con el objetivo de crear correctamente un repositorio para compartir el c´odigo con el director del proyecto y usarlo adecuadamente durante toda la trayectoria del trabajo. La duraci´on de la tarea es de 10 horas. Estudio e implementaci´on de docker (DM2) En esta tarea se realizar´a un estudio del funcionamiento de docker[6] y la implementaci´on de una imagen que contenga todas las dependencias con sus correctas versiones y el c´odigo de la aplicaci´on. Esto nos permite garantizar al usuario que podr´a utilizar correctamente la aplicaci´on. La duraci´on de la tarea es de 20 horas. 17 Revisi´on, an´alisis y ejecuciones del c´odigo en paralelo (DM3) Esta tarea pretende estudiar el ”workflow”, la entrada y salida de datos y las estructuras de datos utilizadas en el c´odigo. Tambi´en se observan los tiempos de ejecuci´on y el consumo de memoria del mismo para la simulaci´on en casos relevantes. Tanto los datos como el estudio previo son necesarios recopilarlos en un documento escrito. La duraci´on de esta tarea es de 40 horas. Refactorizaci´on del c´odigo (DM4) Una vez estudiado el c´odigo esta tarea consiste en el estudiar e introducir varias mejoras en el ”workflow”principal de la aplicaci´on. Esta tarea ser´a realizada en varios sprints de la metodolog´ıa ´agil de Scrum. La duraci´on de esta tarea se estima que sea de 135 horas. Estudio y implementaci´on de las optimizaciones (DM5) Una vez finalizado el desarrollo e implementaci´on de las funcionalidades b´asicas del c´odigo, se iniciar´a la fase de optimizaci´on. En esta tarea se estudiar´a las diferentes herramientas para conseguir un c´odigo limpio y ´optimo (Numba,Cython, etc). La duraci´on de esta tarea es de 80 horas. Ejecuciones en paralelo y an´alisis de rendimiento (DM6) Con una versi´on final del c´odigo se realizar´a un an´alisis del rendimiento con varias ejecuciones en paralelo comparando esta versi´on con la versi´on del c´odigo no optimizado. Su duraci´on es de 35 horas. Empaque (DM7) Finalmente, el c´odigo se subir´a al repositorio del proyecto. ´ Este incluir´a la versi´on final del c´odigo, su documentaci´on y algunos casos de uso. La duraci´on de la tarea se estima que sea de 20 horas. 3.2. Recursos A continuaci´on se muestran todos recursos necesarios para el proyecto. Humanos (R1): El director del proyecto (R1d) y programador (R1p). Materiales (R2): El ordenador personal del programador ya que toda tarea ya sea de gesti´on o del desarrollo se realiza desde casa. 18 Software (R3): Todo el software utilizado es gratuito, entre ellos Git/GitHub, Thonny para la edici´on de c´odigo, docker para la gesti´on de dependencias y LaTeX para la composici´on de documentos de texto. Id Tarea Tiempo Dependencias Recursos GP Gesti´on del proyecto 175h - - GP1 Alcance y contexto 25h - R1,R2 GP2 Planificaci´on temporal 15h GP1 R1,R2 GP3 Presupuesto y sostenibilidad 15h GP1 R1,R2 GP4 Memoria 80h - R1p,R2 GP5 Presentaci´on 20h GP4 R1p,R2 GP6 Reuniones 20h - R1,R2 TP Trabajo Previo 25h - - TP1 Preparaci´on del entorno de trabajo 25h - R1p,R2 DM Desarrollo e implementaci´on de mejoras 340h - - DM1 Estudio de sistema de control de versiones 10h TP1 R1p,R2,R3 DM2 Estudio e implementaci´on de docker 20h TP1 R1p,R2,R3 DM3 Revisi´on, an´alisis y ejecuciones del c´odigo en paralelo 40h TP1,DM1 R1p,R2,R3 DM4 Refactorizaci´on del c´odigo 135h DM3 R1p,R2,R3 DM5 Estudio e implementaci´on de las optimizaciones 80h DM4 R1p,R2,R3 DM6 Ejecuciones en paralelo y an´alisis de rendimiento 35h DM5 R1p,R2,R3 DM7 Empaque 20h DM6 R1p,R2,R3 - Total 540h - - Cuadro 1: Tabla de tareas de la planificaci´on inicial con su duraci´on, dependencias y recursos. Elaboraci´on propia. 19 3.3. Representaci´on gr´afica de la planificaci´on Figura 1: Diagrama de Gantt creado con Gantter. Elaboraci´on propia. 20 3.4. Gesti´on del riesgo Es muy importante preveer los riesgos que se pueden dar durante el proyecto as´ı como contar con alternativas para garantizar que el trabajo finaliza dentro del plazo. Por esa raz´on a continuaci´on se nombran posibles riesgos con sus alternativas. 1. Peque˜nos errores - Estos se tienen en cuenta durante el desarrollo gracias a la metodolog´ıa Scrum ya que nos permitir´a rectificarlos inmediatamente despu´es de detectarlos. Se estima que estos errores consuman 10 horas de trabajo adicionales. 2. Recursos limitados - Como ya se coment´o anteriormente, es posible que las librer´ıas utilizadas por la aplicaci´on nos limiten. Para ello se necesitar´a buscar otras librer´ıas alternativas o incluso desarrollar la funcionalidad necesaria de propia mano. Por ello se cuenta con 1 m´es extra de margen antes de la entrega para poder rectificar pero se estima que esta nueva tarea dure entre 15 y 20 horas. 3. Aver´ıas de recursos materiales - Debido a que el ordenador personal del programador tiene 6 a˜nos y su vida ´util ronda los 8 a˜nos es bastante probable de que sufra una aver´ıa. Si se da el caso, el programador cuenta con un ordenador port´atil que puede sustituir al ordenador principal. Debido a que el ordenador port´atil es mucho menos potente que el ordenador personal puede relentizar el trabajo. Por este motivo contaremos con 15 horas m´as. 3.5. Cambios respecto a la planificaci´on inicial En el momento en que esto es escrito, el proyecto se encuentra en la tarea DM6 - Ejecuciones en paralelo y an´alisis de rendimiento. Por motivos de falta de tiempo y de dificultad a la hora de implementar ciertas funcionabilidades han habido los siguientes cambios respecto a la planificaci´on inicial. 3.5.1. Tareas La decisi´on de volver el repositorio p´ublico de la tarea DM7 - Empaque pasa a estar fuera del trabajo y ser´a decisi´on del Geosciences Applications Group 21 del BSC. Por ello el subobjetivo 6.Crear de un repositorio abierto de f´acil descarga y uso pasar´ıa a ser 6.Crear de un repositorio de f´acil descarga y uso con posiblidad de convertirse en abierto. Id Tarea Tiempo Dependencias Recursos GP Gesti´on del proyecto 175h - - GP1 Alcance y contexto 25h - R1,R2 GP2 Planificaci´on temporal 15h GP1 R1,R2 GP3 Presupuesto y sostenibilidad 15h GP1 R1,R2 GP4 Memoria 80h - R1p,R2 GP5 Presentaci´on 20h GP4 R1p,R2 GP6 Reuniones 20h - R1,R2 TP Trabajo Previo 17h - - TP1 Preparaci´on del entorno de trabajo 17h - R1p,R2 DM Desarrollo e implementaci´on de mejoras 378h - - DM1 Estudio de sistema de control de versiones 10h TP1 R1p,R2,R3 DM2 Estudio e implementaci´on de docker 23h TP1 R1p,R2,R3 DM3 Revisi´on, an´alisis y ejecuciones del c´odigo en paralelo 50h TP1,DM1 R1p,R2,R3 DM4 Refactorizaci´on del c´odigo 155h DM3 R1p,R2,R3 DM5 Estudio e implementaci´on de las optimizaciones 85h DM4 R1p,R2,R3 DM6 Ejecuciones en paralelo y an´alisis de rendimiento 40h DM5 R1p,R2,R3 DM7 Empaque 15h DM1 R1p,R2,R3 - Total 570h - - Cuadro 2: Tabla de tareas de la planificaci´on final con su duraci´on, dependencias y recursos. Elaboraci´on propia. Por lo que hace al resto de tareas se mantienen pero han sufrido cambios en la duraci´on de algunas y en los plazos de tiempo mostrados en el Cuadro 2. Finalmente la duraci´on del trabajo es por lo tanto de 570h. Esto afectar´a en los costes del proyecto explicados m´as detalladamente en la secci´on 6.3.Cambios respecto al presupuesto inicial. En la siguiente secci´on se muestra el diagrama de Gantt de la planificaci´on actual. 22 3.5.2. Representaci´on gr´afica de la planificaci´on actual Figura 2: Diagrama de Gantt creado con Gantter. Elaboraci´on propia. 23 4. Presupuesto 4.1. Identificaci´on de los costes Los costes del proyecto se pueden clasificar en los siguientes tipos: Personales: Director y programador del proyecto. Directos: Ordenador personal y entorno de los trabajadores. Indirectos: Consumo el´ectrico del personal y coste del acceso a internet. 4.1.1. Personales Para calcular los costes de cada persona implicada en el proyecto se cuenta con “Convenio colectivo estatal de empresas de consultor´ıa y estudios de mercado y de la opini´on p´ublica”[8]. Teniendo en cuenta que la jornada m´axima laboral es de 1800 horas anuales podemos realizar los c´alculos representados en la siguente tabla. Perfil Grupo Nivel Salario bru- to(e) SS(e) Coste/a˜no Coste/h Director A - 26790,31 8037,09 34827,40 19,35 Programador C 1 24640,37 7392,11 32032,48 17,80 Cuadro 3: Tabla de los c´alculos de los costes personales. Elaboraci´on propia. Una vez realizados estos c´alculos es posible utilizarlos para calcular los Costes Por Actividad (CPA). 24 “La propiedad industrial e intelectual de los TFG de modalidad A est´a regulada por la normativa aprovada por el Consejo de Gobierno(10/10/2008) por la qual se aprueba la confidencialidad, responsabilidad patrimonial y propiedad indistrial e intelectual a la UPC. Corresponder´a a la UPC la titularidad sobre las invenciones desarrolladas exclusivamente por los estudiantes si se ha desarrollado en el marco de una actividad acad´emica que ha estado dirigida y/o coordinada por el profesorado de la UPC. En el caso que el desarrollo de la obra intelectual haya estado dirigida y/o coordinada por el profesorado de la UPC, corresponder´a a la UPC la titularidad de los derechos de explotaci´on sobre esta obra y el estudiante y el profesor ser´an considerados coautores de la misma. En el caso de explotaci´on de la obra por parte de la UPC que le suponga un beneficio econ´omico, el autor o conjunto de autores tendr´an derecho a una participaci´on del 50 % de los beneficios netos obtenidos.” 6.2. Licencias Todas las licencias del software utilizado para el desarrollo de este proyecto son de c´odigo abierto y libre. 31 7. Propagador de ondas ac´usticas 2D/3D (PSM) La precisi´on y eficiencia de las t´ecnicas de im´agenes geof´ısicas como inver- si´on de forma de onda completa (FWI) o migraci´on de tiempo inverso (RTM) dependen en gran medida de las t´ecnicas num´ericas utilizadas para modelar la propagaci´on de ondas s´ısmicas. En la aproximaci´on ac´ustica, las fases se pueden reproducir correctamente mediante la ecuaci´on de onda escalar. Existen muchos m´etodos para modelar la propagaci´on de ondas en medios ac´usticos. Los m´etodos m´as comunes son el expl´ıcito en el tiempo y el orden superior en diferencias finitas espaciales (FD). A pesar del amplio ´exito de este m´etodo, la aproximaci´on precisa y eficiente de las derivadas de tiempo y espacio continua siendo un tema abierto al estudio. Por un lado, tradicionalmente la integraci´on de tiempo se obtiene utilizando esquemas de FD de bajo orden. Por otro lado, la precisi´on espacial podr´ıa mejorarse mediante el uso de diferenciaci´on espectral. El m´etodo de dominio de tiempo pseudoespectral de Fourier (PSTD) se basa en dicha diferenciaci´on espectral por medio de la transformada de Fourier (TF). Dadas sus ventajas en precisi´on num´erica, el c´odigo utilizado y mejorado en este TFG se basa en el m´etodo PSTD. El problema f´ısico considerado en este trabajo est´a governado por la ecuaci´on ac´ustica ∂2u(x, t) ∂t2+L2u(x, t) = 0 (1) donde u(x, t) es la presi´on ac´ustica en el tiempo t. El operador espacial se define de acuerdo a L=lc∇donde l=√−1. Se considera un dominio computacional x∈Ω en 2D y 3D a trav´es del cual la velocidad del sonido c(x) podr´ıa variar. Para el t´ermino espacial continuo de la ecuaci´on 1, el t´ermino L2puede expandirse como L2u(x, t) = −c2(x){F−1 x[k2 xFx[u(x, t)]]+F−1 y[k2 yFy[u(x, t)]]+F−1 z[k2 zFz[u(x, t)]]} (2) 32 donde kx,ky, y kzson el n´umero de onda en cada dimensi´on cartesiana; FyF−1denotan la TF continua en 1D y su inversa, respectivamente. La formulaci´on de los m´etodos de Fourier PS se basa en un an´alogo discreto de la ecuaci´on 2, en el que las derivadas 1D se calculan en el dominio del n´umero de onda mediante TF discreto (DFT) y se devuelven al dominio espacial mediante una TF inverso. Por ello, se ha considerado un medio rectangular Ω discretizado por una cuadr´ıcula regular de igual espaciado ∂, es decir, x=i∂,y=j∂,z=l∂, y denotamos los puntos de la cuadr´ıcula por x. Adem´as, denotamos como L2la aproximaci´on PS de L2en un muestreo de cuadr´ıcula, de acuerdo con la ecuaci´on 2. Por lo tanto, L2es un operador de matriz discreta con las mismas dimensiones que la malla rectangular. Es importante se˜nalar que la DFT es una buena aproximaci´on a la FT siempre que se cumpla el teorema de muestreo de Nyquist-Shannon en cada dimen- si´on, por ejemplo, kmax ≤π/∂, y la extensi´on peri´odica de la distribuci´on espacial es continua. En este caso, el esquema PS tiene una precisi´on espacial de O(∂N). 7.1. Algoritmo de migraci´on de tiempo inverso El m´etodo num´erico descrito anteriormente se ha utilizado para implementar el algoritmo de migraci´on de tiempo inversi´on (RTM). RTM es un una migraci´on de ecuaci´on de onda ac´ustica bidireccional preapilada para obtener im´agenes precisas en y debajo de ´areas con grandes complejidades estructurales y de velocidad constante/variable. RTM tiene un historial probado para generar im´agenes de alta calidad y se utiliza cada vez m´as para refinar los l´ımites estructurales de dominios computacionales durante la construcci´on de modelos de velocidad o de presi´on. Para obtener im´agenes de alta fidelidad, RTM calcula soluciones num´ericas para la ecuaci´on de onda completa (ver ecuaci´on 1). Como tal, no tiene limitaci´on de inmersi´on y maneja todas las formas de onda complejas de m´ultiples rutas. Hist´oricamente, RTM se consider´o poco pr´actico debido a los altos costos computacionales y una mayor sensibilidad a los par´ametros de velocidad y reflectividad que los m´etodos continuos unidireccionales m´as establecidos. Sin embargo, la mejora constante de las arquitecturas computaciones, junto con flujos de trabajo de creaci´on de modelos de velocidad/presi´on m´as precisos, han hecho viable el desarrollo e implementaci´on de algoritmos RTM para la generaci´on de im´agenes en diferentes medios. 33 7.2. Estructura del c´odigo El prototipo resuelve el problema de la propagaci´on de una onda ac´ustica a trav´es de diferentes medios. Se ejecuta con un peque˜no script llamado “run” que se encarga de lanzar tres ejecuciones, una para cada fase del programa. Estas tres fases se denominan: preprocesado,kernel ypostprocesado. 7.2.1. Dependencias Primeramente para poder ejecutar el prototipo es necesaria la preparaci´on del entorno de desarrollo con la instalaci´on de dependencias y software de terceros. A continuaci´on se expondr´an todas las herramientas utilizadas. El sistema operativo utilizado es linux en concreto he escogido la distribuci´on de Ubuntu en su versi´on 18.04 LTD. Una vez instalado el sistema operativo se requieren instalar las dependencias. MPICH es una implementaci´on portable, de alto rendimiento y open source de MPI, un est´andar para el env´ıo de mensajes para aplicaciones con memoria distribuida utilizada en la computaci´on paralela. HDF5 es un formato de almacenamiento de datos jer´arquico pensado para organizar grandes cantidades de datos. Tiene integrado herramientas de al- to rendimiento para optimizar tanto el tiempo de acceso a los datos como el espacio de almacenamiento. Adem´as, tambi´en ofrece herramientas para consultar, manipular, visualizar y analizar los datos. FFTW es conocida por ser la implementaci´on open source m´as rapida de la “Fast Fourier Transform (FFT)”, un algoritmo eficiente que permite calcular la transformada de fourier discreta (DFT) en una o m´as dimensiones de tama˜no de entrada arbitrarios. A continuaci´on se instal´o Python en su versi´on 3.8.5, pero al tener problemas a la hora de instalar ciertos paquetes de Python necesarios finalmente se opt´o por la versi´on 3.6.9. Los paquetes de Python utilizados son los siguientes. Numpy: proporciona objetos de arrays multidimensionales conjunatamente con una amplia lista de rutinas muy eficientes para operar dichos arrays. 34 mpi4py: nos permite utilizar MPI desde Python. mpi4py-fft: incluye una interfaz de FFTW para Python y junto a mpi4py nos permite distribuir arrays muy grandes gesationando las comunicaciones entre procesos h5py: es una interfaz de HDF5 para Python que nos permite almacenar grande cantidades de datos num´ericos y facilidades para manipular los datos desde Numpy. PyYAML: YAML es un formato de serializaci´on de datos dise˜nado para la legibilidad humana. PyYAML es un parser y emisor para Python de YAML. Cython: Es un lenguaje de programaci´on para simplificar la escritura de m´odulos de extensi´on para Python en C y C++. Scipy: Es una biblioteca libre y de c´odigo abierto que se compone de herramientas y algoritmo matem´aticos. Sphinx: Librer´ıa para generar documentaci´on de Python. singleton decorator: Permite crear objetos “singleton” a partir de un decorador. colorama: Permite colorear el texto en la terminal y el posicionamiento del cursor. matplotlib: Biblioteca para visualizaci´on gr´afica en 2D y 3D. Al momento de tener todas las dependencias necesarias para ejecutar el prototipo se estudi´o el funcionamiento de Git y GitHub para poder crear un repositorio en el que se pudiese compartir el prototipo y as´ı en un futuro ir actualiz´andolo con las mejoras implementadas. 7.2.2. Entrada y salida La entrada se compone del archivo “params.yaml” que contiene los par´ametros mostrados en la Figura 3. 35 Figura 3: params.yaml Adem´as hay que tener un archivo HDF5 con los datos del modelo de velocidades, la fuente y el receptor. Por ello contamos con dos scripts escritos en Python “build IR 2Ddata.py” y“build IR 3Ddata.py” que se encargan de crear los archivos “IR input data2D.h5” y“IR input data3D.h5” respectivamente. Estos contienen un modelo de velocidad homog´eneo, es decir, la velocidad es constante para todos los puntos del modelo y por lo tanto afectar´an de la misma forma a la onda. La salida de la aplicaci´on es un archivo HDF5 con los siguientes datos: axis order: Orden de los ejes (ZX/ZXY). code version: Versi´on del c´odigo. date: La fecha de la ejecuci´on. machine: Nombre de la m´aquina en la que se realiz´o la ejecuci´on. dh: Variaci´on en el espacio (Delta in space). 36 dt: Variaci´on en el tiempo (Delta in time). nt: N´umero de pasos temporales. num dimensions: N´umero de dimensiones (2D/3D). num grid points: N´umero de puntos de la cuadr´ıcula. num processors: N´umero de procesos con los que se ha ejecutado el programa. image: Valores de la imagen de la propagaci´on de la onda captada por el receptor. 7.2.3. Preprocesado El preprocesado es totalmente secuencial y est´a compuesto de las siguientes partes. 1. Lectura de la entrada: Consiste en parsear el archivo de entrada y crear objetos con todos los datos del archivo. 2. Inicializaci´on: Esta parte se encarga inicializar el conjunto de datos del modelo, las fuentes y los receptores haciendo lecturas del archivo hdf5 especificado en la entrada. Adem´as, se crea el directorio temporal y el directorio donde se alojar´a el archivo de salida. 3. Calcular la cuadr´ıcula espacio-tiempo: En esta parte se calcula la cuadr´ıcula espacio-temporal que es una discretizaci´on del modelo de entrada y por lo tanto se interpolar´a el modelo a dicha cuadr´ıcula. Esta discretizaci´on est´a adaptada al proceso, es decir que el tama˜no de la cuadr´ıcula depender´a del modelo, la frecuencia m´axima de la onda y la variable de entrada de ppw (points per wavelength). En la Figura 4 se muestra un ejemplo. Se guardan tanto los datos como los metadatos calculados en archivos temporales. 4. Calcular metadatos del receptor: Consiste en computar los datos necesarios para el receptor y en guardarlos en archivos temporales. 37 Figura 4: Ejemplo de cuadr´ıcula espacio-tiempo para un modelo 2D (izquierda) y boceto de cuadr´ıcula con ppw=2 (derecha). Elaboraci´on propia 7.2.4. Kernel El kernel es la ´unica parte del prototipo que se hace en paralelo y se compone de las siguentes partes. 1. Lectura de la entrada: Es exactamente la misma parte que en el preprocesado. 2. Inicializaci´on: Lo ´unico que lo diferencia de la fase del preprocesado es que se inicializan los timers. 3. Importar metadatos de la cuadr´ıcula espacio-tiempo: Se cargan los metadatos de la cuadr´ıcula espacio-temporal calculada en la fase del preprocesado leyendo el archivo. 4. Importar metadatos del receptor: Se cargan los datos del receptor calculada en la fase del preprocesado leyendo el archivo. 5. Crear estructuras paralelas: Creaci´on de las estructuras de datos necesarias para paralelizar la aplicaci´on y lectura de los datos de la cuadr´ıcula. 38 6. Crear estructuras de la fuente y del receptor: Creaci´on de las estructuras de la fuente y calcular la descomposici´on del dominio del receptor 7. Crear estructuras de los l´ımites: Consiste en inicializar los datos de los l´ımites. 8. Inicializar y ejecutar solver: Inicializa los datos del Solver y lo ejecuta. El Solver es el que se encarga de resolver el sistema de ecuaciones y de toda la parte matem´atica del problema. Cuenta con dos fases, “forward” donde se simula hacia delante en el tiempo (0..nt) y “backward” donde se simula hacia atr´as en el tiempo (nt..0). 7.2.5. Postprocesado Por ´ultimo la fase del postprocesado que tambi´en es completamente secuencial y est´a compuesto de las siguientes partes. 1. Lectura de la entrada: Es exactamente la misma parte que en el preprocesado. 2. Inicializaci´on: Es exactamente la misma parte que en el preprocesado. 3. Importar metadatos de la cuadr´ıcula espacio-tiempo: Es exactamente la misma parte que en el kernel. 4. Postprocesar Soluci´on: Lectura de los datos calculados en el Solver, interpolaci´on de la cuadr´ıcula al modelo y escritura de los datos de salida. El prototipo cuenta con dos modos de ejecuci´on “modelling” e“image”. El modo “modelling” genera un fichero que puede ser dado como entrada al modo “image” y tan solo ejecuta la fase forward del Solver. Sin embargo en el modo “image” se ejecutan las dos fases. 39 7.3. Simulaciones y resultados 7.3.1. Test de validaci´on: modelo impulso-respuesta Para verificar la implementaci´on del algoritmo RTM, se una utilizado un modelo de impulso-respuesta (I-R). Una funci´on de I-R describe la evoluci´on de la variable de inter´es a lo largo de un horizonte de tiempo especificado despu´es de un choque de onda en un momento dado. En otras palabras, una funci´on I-R define c´omo un sistema/modelo responde a alguna se˜nal de entrada o impulso. Una descripci´on amplia y detallada se puede encontrar en Fletcher, R. (2009)[11]. Este test lo utilizaremos para verificar la precisi´on num´erica y performance (escalabilidad y tiempos). El modelo I-R utilizado cuenta con tres impulsos en tiempo (tres picos). Los detalles de cada configuraci´on y experimento se detallan en secciones m´as adelante. 7.3.2. Precisi´on num´erica Para verificar que los resultados en secuencial son correctos contamos con tres archivos binarios que provisionan de la imagen, la fuente y el modelo de referencia de la salida del caso base homog´eneo 2D. Adem´as contamos con un script para visualizar tanto nuestra salida, como la de referencia y sus diferencias. 40 Simulaciones 3D Figura 10: Observable 1 (3D). Elaboraci´on propia. Por lo que hace a las simulaciones 3D podemos observar los mismos resultados. En la Figura 10 vemos que con el modo “image” las fases que peor escalan son la comunicaci´on y la IO por las razones comentadas anteriormente. 47 Figura 11: Observable 2 (3D). debug snapshot time skip = 2 (izquierda) y debug snapshot time skip = 6 (derecha). Elaboraci´on propia. En la Figura 11, la IO se comporta como se esperaba a la hora de cambiar el valor del par´ametro “debug snapshot time skip”. Figura 12: Observable 3 (3D). Elaboraci´on propia. 48 Sin embargo, en la Figura 12, a diferencia de las simulaciones 2D esta vez las comunicaciones se mantienen practicamente constantes para el modo image. Y para el modo “modelling” la fase otros tiene mucho peor rendimiento para casos grandes como es la simulaci´on 3D. 8. Implementaci´on de mejoras 8.1. Docker Al contar con incompatibilidades con algunas versiones de algunas de las dependencias mencionados anteriormente se tom´o la decisi´on de utilizar docker para remediarlo. Docker permite a los usuarios crear aisladamente e independientemente entornos para ejecutar y lanzar sus aplicaciones. Estos entornos se llaman contenedores. Con esto podemos crear un contenedor con todo el software necesario y con sus versiones compatibles para as´ı poder ejecutar nuestra aplicaci´on en cualquier m´aquina que disponga de docker. Para crear este contenedor es necesario crear una imagen de docker a partir de un archivo creado por nosotros llamado “Dockerfile”. Con este archivo podemos utilizar una imagen creada por docker de Ubuntu 18.04 LTD e instalar todas las dependencias. Una vez hecha la imagen es posible ejecutarla y hacer uso de nuestra aplicaci´on. 8.2. Flujo de trabajo del c´odigo El primer cambio consisti´o en actualizar las funciones de medici´on de tiempos. M´as concretamente, se han implementado timers de MPI en lugar de timers de funciones nativas de Python. Esto permite obtener mediciones de tiempo m´as precisas en ejecuciones paralelas. A partir de observar qu´e hay partes repetidas entre las tres fases hab´ıa que integrarlas, manteniendo as´ı en memoria principal los datos y objetos que necesitamos permiti´endonos eliminar lecturas y escrituras para el paso de datos y metadatos entre las fases. 49 Una vez integradas las fases, el siguiente paso fue paralelizar tanto la fase del preprocesado como la fase del postprocesado para aprovechar al m´aximo el paralelismo. Para ello se necesita paralelizar las interpolaciones de las dos fases distribuyendo as´ı el modelo y la cuadricula entre los cores. Gracias a la librer´ıa de “mpi4py fft” podemos utilizar su m´odulo PFFT para distribuir las matrices de dos formas diferentes, Slab decomposition donde tan solo una dimensi´on est´a distribuida y Pencil decomposition donde dos dimensiones est´an distibuidas. Con esta librer´ıa estamos obligados a tener como m´ınimo un eje alineado lo que implica que para 2D tan solo podemos utilizar la Slab decomposition. En las Figuras 13 y 14 se muestra un ejemplo de cada descomposici´on. Figura 13: Ejemplo de slab decomposition con 4 cores. (a) descomposici´on en el eje Y; (b) descomposici´on en el eje X. Tomado de la p´agina de 2decomp&FFT[12]. Figura 14: Ejemplo de pencil decomposition con 4*3 cores. (a) X-pencil; (b) Y-pencil; (c) Z-pencil. Tomado de la p´agina de 2decomp&FFT[12]. 50 Si solo contamos con modelos de velocidad homog´eneos no se requiere ninguna implementaci´on adicional. Sin embargo, este no es el caso del problema que se pretende resolver en este TFG. El programa debe aceptar escenarios realistas y complejos, por lo que supone que los modelos de velocidad ser´an en su mayor´ıa heterog´eneos. Para conseguir una interpolaci´on en paralelo correcta y suave de un modelo de dichas caracter´ısticas, es necesario que cada proceso se comunique con sus vecinos para obtener la serie de datos que requieren y con ellos poder computar correctamente los valores de la frontera de cada proceso. Para esto se a˜nadi´o un par´ametro interpolation: halo que indica la anchura de puntos consultados a cada vecino. Ejemplo en la Figura 15. Figura 15: Ejemplo de comunicaci´on para la interpolaci´on con halo = 3. Elaboraci´on propia. Para comprobar que la implementaci´on de la interpolaci´on es correcta se han implementado test unitarios que nos permiten cambiar el tama˜no de la entrada y de la salida adem´as de sus valores homog´eneos o heterog´eneos(randoms). Pero esta implementaci´on nos acaba dando problemas a la hora de utilizarla con muchos cores y un gran consumo de memoria, ya que las funciones para interpolar que utilizamos de la librer´ıa Scipy consumen bastante memoria. Por lo que acabamos implementando una alternativa pero esta vez en secuencial. 51 La idea es dividir la matriz de valores a interpolar en partes peque˜nas teniendo en cuenta que necesitamos el halo de valores de las partes vecinas y iterar sobre ellas realizando interpolaciones peque˜nas y acabar haciendo una reducci´on para construir el resultado. El tama˜no de cada iteraci´on se ve definido por otro par´ametro introducido llamado interpolation: slice width. 9. Conclusiones y trabajo a futuro En este TFG se han estudiado e implementado mejoras a un prototipo PSM para la simulaci´on de ondas ac´usticas en 2D/3D. Se han analizando y desarrollado diferentes estrategias computacionales impactando positivamente en la eficiencia del m´etodo num´erico y en la flexibilidad del c´odigo. Todo lo anterior bajo un enfoque open source. La conclusi´on principal es que los objetivos se han alcanzado satisfactoriamente. Se ha demostrado que para cada una de las fases del algoritmo se ofrece un rendimiento distinto y se han implementado mejoras sustanciales en el c´odigo. Adem´as, se ha creado un repositorio con soporte para Dockers que resultar´a de utilidad para futuros desarrollos a partir de la versi´on implementada en este TFG. A nivel personal me ha aportado nuevos conocimientos t´ecnicos y a su vez me ha permitido aportar y aplicar al proyecto todo el conocimiento adquirido durante los a˜nos de este grado. Tambi´en me ha aportado nuevos conocimientos para poder desarrollar exitosamente toda la parte de gesti´on de un proyecto. No deja de ser cierto que, aunque la experiencia del proyecto ha sido enriquecedora y que la pandemia del SARS-CoV 2 no lo ha permitido, me hubiese gustado vivir la experiencia de la forma usual. Como hemos comentado anteriormente, el trabajo a futuro consiste en mejorar la tarea de IO del algoritmo ya que el algoritmo FWI presentado en este documento requiere escrituras a disco eficientes. Para ello se podr´ıa implementar y analizar escrituras con un formato diferente como podria ser en binario, o buscar m´etodos para escribir menos bytes e incluso utilizar un sistema de compresi´on de datos. Tambi´en se podr´ıan utilizar nuevas tecnolog´ıas de memoria como Optane o aprovechar escrituras en memoria virtual. 52 A parte de esto, ser´ıa bueno extrapolar el tiempo por iteraci´on para aproximar el tiempo total con varias mallas de tama˜no representativo. Tambi´en cambiar el repositorio de privado a p´ublico una vez se crea oportuno que est´e listo para ser open source. 53 Referencias [1] Geosciences Applications Group of the Barcelona Supercomputing Center. https://www.bsc.es/es/discover-bsc/organisation/ scientific-structure/geophysical-applications [2] Octavio Castillo Reyes. https://futur.upc.edu/OctavioCastilloReyes [3] Google Meet. https://meet.google.com/ [4] Git. https://git-scm.com/ [5] GitHub. https://github.com/ [6] Docker. https://www.docker.com/ [7] Ubuntu. https://ubuntu.com/ [8] Convenio colectivo estatal de empresas de consultor´ıa y estudios de mercado y de la opini´on p´ublica. https://www.boe.es/boe/dias/2018/03/06/pdfs/ BOE-A-2018-3156.pdf [9] Tabla de coeficientes de amortizaci´on lineal. https://www.agenciatributaria.es/AEAT.internet/Inicio/ _Segmentos_/Empresas_y_profesionales/Empresas/Impuesto_ sobre_Sociedades/Periodos_impositivos_a_partir_de_1_1_ 2015/Base_imponible/Amortizacion/Tabla_de_coeficientes_ de_amortizacion_lineal_.shtml [10] Normativa TFG GEI FIB UPC. https://www.fib.upc.edu/sites/fib/files/documents/estudis/ normativa-tfg-mencio-addicional-gei-br.pdf 54 [11] Fletcher, R. P., Du, X., & Fowler, P. J. (2009). Reverse time migration in tilted transversely isotropic (TTI) media. Geophysics, 74(6), WCA179- WCA187. [12] Library for 2D pencil decomposition and distributed Fast Fourier Transform. http://www.2decomp.org/decomp.html 55 Anexo Anexo A: Test Escalabilidad Este anexo contiene todas las gr´aficas obtenidas a partir de pruebas de escalabilidad ejecutadas en MN. Consultalas en las siguientes p´aginas ↓↓↓. 56 Figura 22: Gr´aficas del tiempo de ejecuci´on (izquierda) y speedup (derecha) del Solver 3D con el modo “modelling” con debug snapshot time skip = 2. Elaboraci´on propia. 63 Figura 23: Gr´aficas del tiempo de ejecuci´on (izquierda) y speedup (derecha) del Solver 3D con el modo “modelling” con debug snapshot time skip = 6. Elaboraci´on propia. 64