Modelado Sísmico Computacional: resolución mediante elementos finitos en entornos de computación en paralelo

Los objetivos generales que nos hemos impuesto en esta memoria son los siguientes:

1. Formular completamente el problema de la elastodinámica lineal, enunciando todas las hipótesis y principios fundamentales que, desde la Mecánica de los Medios Continuos, permitan un acercamiento global al modelado de ondas elásticas. Desde este punto de vista, deberán describirse las ecuaciones en derivadas parciales (PDEs) que describen el movimiento de un medio heterogéneo.

2. Deducir y estudiar la formulación variacional que permite discretizar mediante elementos finitos el problema de la propagación de ondas elásticas. Para ello habremos hacer uso de las herramientas que provee el análisis funcional, y se enunciarán las condiciones que justifican la unicidad y existencia de la solución. En este punto, los métodos del espacio de Sobolev habrán de desempeñar un papel crucial en su análisis.

3. Aplicar el método de los elementos finitos al problema de la elastodinámica formulado variacionalmente. Exponer el error derivado de la discretización y fundamentar las ventajas e inconvenientes del método. A su vez, estudiar y encontrar las técnicas adecuadas para los problemas asociados al mismo: representación de la fuente, condiciones de amortiguación en los límites del dominio, diagonalidad de las matrices, etc.

4. Encontrar los métodos numéricos más adecuados que integren la ecuación diferencial ordinaria (formulada matricialmente y con coeficientes constantes) que se obtiene al discretizar el problema con elementos finitos. Elegir las técnicas más convenientes para resolver el sistema de ecuaciones final. Fundamentar los criterios de elección de esquemas explícitos o implícitos.

5. Estudiar con detalle algún problema sismológico concreto. Realizar simulaciones en modelos que se hayan resuelto anteriormente con otros métodos de discretización, poniendo particular atención en aquellos que puedan crear dificultades. Debido a la ausencia de datos reales provenientes del campo de la exploración sísmica con los que comparar los resultados, repetir en lo posible simulaciones sobre medios complejos que se puedan hallar en la literatura.

6. Comprender las necesidades y características de la computación en sistemas multiprocesadores. Estudiar las diferentes arquitecturas de ordenador, los paradigmas más importantes de programación en paralelo, y las métricas adecuadas que permitan establecer consecuencias y conclusiones de su uso.

7. Estudiar los métodos multidominio que den lugar a una implementación paralela eficiente del problema de propagación de ondas. Analizar en detalle los aspectos teóricos y matemáticos de las técnicas de descomposición del dominio. Proponer y enunciar un método que mejore las ventajas de los actuales métodos de subestructuración.

8. Realizar simulaciones en paralelo en sistemas sencillos mediante las técnicas estudiadas en el apartado anterior. Deberá hacerse especial hincapié en las diferencias que pudieran surgir entre los algoritmos secuencial y paralelo, tomando el primero como medida relativa del error de la aproximación. Simular campos de ondas en medios geológicos muy grandes que requieran enormes tiempos de computación. Comparar la técnica propuesta con las ya existentes.

En conexión con estos propósitos, la presente memoria se estructura de la siguiente manera:

• Un primer capítulo de introducción general al modelado sísmico computacional, antecedentes y descripción del estado actual del uso de supercomputación para la resolución del problema de la propagación de ondas sísmicas.

• Tres capítulos en los que se desarrolla ampliamente toda la teoría: formulación del problema de la elastodinámica lineal a partir de los principios fundamentales de la Mecánica de los Medios Continuos y las leyes del comportamiento de los materiales; formulación variacional del problema y resultados de existencia y unicidad de la solución; aplicación del método de los elementos finitos al problema de propagación de ondas elásticas. Formulación matricial. Ventajas e inconvenientes. Condiciones de contorno absorbentes y descripción de la fuente sísmica empleada.

• Un quinto capítulo en el que se proponen métodos numéricos eficientes para la integración temporal de la ecuación diferencial de segundo grado, estudio de la exactitud de dichos esquemas, métodos de resolución de ecuaciones y técnicas de almacenamiento de matrices sparse.

• Un capítulo de resultados con diversas simulaciones dé campos de ondas a través de diferentes medios heterogéneos.

• Dos capítulos dedicados a introducir los conceptos esenciales de la computación en paralelo, las diferentes arquitecturas y los métodos de descomposición del dominio cuyo diseño responde a las características buscadas para realizar esta memoria.

• Un último capítulo dedicado a las simulaciones efectuadas mediante paralelismo y métodos de descomposición del dominio.

Debido a las abundantes referencias que surgen a lo largo de la presente memoria, la bibliografía aparece anexa a cada capítulo.

En el principio… 5

Objetivos y descripción de la memoria 7

1 Introducción general 11
   1.1 El modelado sísmico 11
         1.1.1 Antecedentes no numéricos 12
         1.1.2 Métodos numéricos 13
         1.1.3 Tres ejemplos de discretización 15
   1.2 Antecedentes 19
         1.2.1 Métodos interiores 19
         1.2.2 Métodos de contorno 24
         1.2.3 Métodos híbridos 26
Bibliografía 26

I Propagación de ondas mediante elementos finitos 33

2 Formulación del problema 35
   2.1 Hipótesis 36
   2.2 Principios fundamentales 36
         2.2.1 Principio de conservación de la masa 36
         2.2.2 Principio fundamental de la .dinámica 37
         2.2.3 Término de fuente 38
    2.3 Leyes del comportamiento 38
        2.3.1 Elasticidad lineal isótropa 39
        2.3.2 Anisotropía 39
        2.3.3 Conservación de la energía 40
        2.3.4 Atenuación 41
   2.4 Condiciones iniciales y de contorno 42
        2.4.1 Equivalencia entre los problemas de Cauchy y de Green 43
Bibliografía 44

3 Formulación variacional 47
   3.1 Espacios de Sobolev 47
        3.1.1 Normas 48
        3.1.2 Tres resultados de gran utilidad 49
   3.2 Ecuación fundamental de la elastodinámica 49
        3.2.1 Existencia y unicidad de la solución 52
Bibliografía 54

4 El Método de los Elementos Finitos 55
   4.1 El problema discreto 55
        4.1.1 Estudio del error de la aproximación 57
   4.2 El problema matricial 58
        4.2.1 Concentración o diagonalización de matrices 60
   4.3 Ventajas e inconvenientes 62
   4.4 Condiciones de contorno absorbentes 64
       4.4.1 Zonas absorbentes 66
       4.4.2 Amortiguador lineal modificado 69
   4.5 La fuente sísmica 72
       4.5.1 Función temporal 73
       4.5.2 Función espacial 75
Bibliografía 76

5 Métodos numéricos 81
   5.1 Métodos de integración en el tiempo 82
         5.1.1 Métodos de Newmark 82
         5.1.2 Método de Houbolt 86
         5.1.3 Método de Wilson-9 86
   5.2 Exactitud de los esquemas de integración 87
   5.3 Métodos de resolución de ecuaciones 88
         5.3.1 Métodos directos 88
        5.3.2 Métodos iterativos 89
        5.3.3 ITPACKV-2D 90
   5.4 Matrices sparse 91
        5.4.1 Compressed Sparse Row 92
        5.4.2 Formato diagonal 92
        5.4.3 Formato ELLPACK-ITPACK 93
        5.4.4 Formato skyline 93
Bibliografía 94

6 Resultados 97
   6.1 Interfases entre medios isótropos 98
         6.1.1 La interfase líquido-sólido 98
         6.1.2 Interfases sólido-sólido 99
   6.2 Modelo de inclusión de hidrocarburos 103
   6.3 Modelo de intrusión 109
   6.4 El modelo cuña 112
   6.5 El domo salino 113
   6.6 El modelo Marmousi 121
Bibliografía 130

II Paralelización del problema de propagación de ondas elásticas 133

7 Paralelismo como solución 135
   7.1 Estructura de un método DDM paralelo 137
   7.2 La computación en paralelo 138
         7.2.1 Ley de Amdahl 140
         7.2.2 Parámetros importantes 141
         7.2.3 Equilibrio en la ejecución 142
   7.3 Introducción a las arquitecturas paralelas 142
         7.3.1 Single Instruction Multiple Data 143
         7.3.2 Multiple Instruction Multiple Data 143
   7.4 Conceptos de programación paralela 147
         7.4.1 Message-Passing 147
         7.4.2 Data Parallel 150
Bibliografía 151

8 Técnicas de descomposición 153
   8.1 Métodos de descomposición del dominio 153
   8.2 El método de Schwarz 154
   8.3 El método del complemento de Schur 157
   8.4 Método multidominio de Quarteroni 159
   8.5 Análisis 163
        8.5.1 DDM y propagación de ondas 163
        8.5.2 El complemento de Schur 164
        8.5.3 El problema de los sistemas algebraicos explícitos 165
   8.6 El método multibloque 166
        8.6.1 Ventajas 170
        8.6.2 Desventajas 171
Bibliografía 172

9 Resultados 175
   9.1 Entorno data-parallel 175
         9.1.1 Análisis de eficiencias en data-parallel 176
         9.1.2 Solución PowerFortran 177
   9.2 Message-passing 179
   9.3 Esquemas explícitos paralelos 185
        9.3.1 Método de subestructuración 185
        9.3.2 Método multibloque 187
        9.3.3 Comparación y conclusiones 190
   9.4 Esquemas implícitos paralelos 191
        9.4.1 Método de iteración por subdominios 194
        9.4.2 Método multibloque 196
        9.4.3 Comparación y conclusiones 207
Bibliografía 209

Conclusiones 212

2. Se ha aplicado el método de los elementos finitos a la formulación variacional empleando la semidiscretización de Galerkin. El problema se ha planteado en su generalidad, incluyendo las adecuadas condiciones absorbentes y el término de fuente que dotan de sentido físico a la solución.
A este respecto, la elección de una zona viscoelástica alrededor del dominio de interés, con un coeficiente de amortiguamiento que varía linealmente en el dominio de definición, se ha percibido como la más sencilla y efectiva.
De las diferentes posibilidades de fuentes sísmicas consultadas, todas producen buenos resultados al modular correctamente las frecuencias de interés, si bien en nuestro trabajo hemos optado por definir la fuente también en el dominio espacial de manera que los fenómenos de dispersión de frecuencias no admisibles por la malla sean mínimos.

3. Se ha analizado las posibilidades de emplear esquemas explícitos o implícitos para la resolución del problema. En resumen, se trata de comprobar si la diagonalización de las matrices que intervienen en el problema es o no conveniente. En este punto, prácticamente no hay unanimidad en la literatura consultada, si bien en nuestras pruebas el método implícito de aceleración promedio constante ha proporcionado resultados más satisfactorios. El esquema de Diferencias Centrales con matrices de masa y amortiguamiento construidas diagonalmente, dan lugar a simulaciones baratas y rápidas, pero su exactitud no es siempre suficiente al tener que determinar con exactitud el rango de estabilidad del esquema.
Se ha empleado un resolutor de ecuaciones de contrastada calidad y optimización, el ITPACKV-2D. Los objetivos perseguidos en esta memoria no suponían la programación de métodos del subespacio de Krylov, habida cuenta de su abundancia tanto en la literatura especializada como en Internet. Se ha desechado el empleo de resolutores directos por su menor eficiencia con grandes matrices sparse y mayor costo de almacenamiento.

4. Se han resuelto diversos problemas de propagación de ondas, especialmente aquellos que dan lugar a dificultades en otros esquemas de discretización. La ausencia de datos empíricos de exploración de hidrocarburos nos ha obligado a ejecutar un proceso de propagación de ondas sísmicas en aquellos modelos que aparecen en las referencias bibliográficas más importantes. Asimismo, se ha procedido a simular un frente de ondas que avanza por el modelo de Marmousi como demostración de que el método de los elementos finitos produce tiempos de propagación satisfactorios y contrastados.

5. Se han estudiado las técnicas de descomposición del dominio para grandes estructuras que pudieran ser más efectivas en el caso de propagación de ondas sísmicas. Este tipo de técnicas son aún poco conocidas en la comunidad de sismólogos con interés en el modelado, pero se trata de herramientas muy potentes y adecuadas para el cálculo en entornos de computación en paralelo. Se ha demostrado que las técnicas multidominio producen eficiencias mejores a las de descomposición de los datos al aplicarse al método de los elementos finitos.

6. Se ha enunciado y formulado una técnica de descomposición del dominio nueva y original basada en la equivalencia de los problemas de Cauchy y de Green para la elastodinámica lineal. Esta técnica, denominada multibloque, se ha formulado atendiendo las necesidades de conseguir un cálculo rápido en esquemas implícitos, pues los métodos habituales (complemento de Schur, método multidominio) no resultan ni muy eficientes ni sencillos de programar.

7. Se ha comparado la técnica multibloque con la de subestructuración para el caso de esquemas explícitos. Se ha comprobado que en sistemas de memoria distribuida, resulta aquella más eficiente cuando el número de procesadores es alto. En entornos de memoria compartida, la técnica multibloque es siempre más eficiente a igualdad de procesadores. Los errores que se cometen con ambas técnicas son mínimos (respecto al problema secuencial).

8. Se ha implementado un esquema de iteración por subdominios para intentar resolver el problema implícito. El método no ha convergido.

9. Se ha implementado con éxito la técnica multibloque en esquemas implícitos, habiéndose encontrado una eficiencia muy alta en los resultados obtenidos. Los errores que producen son admisibles habida cuenta que, el acoplamiento entre frontera e interior de un subdominio, da lugar a la aparición de reflexiones espúreas de baja amplitud, en torno al 7% respecto la fase principal, posiblemente achacable a la precisión del cálculo de la máquina.

10. Para la realización del trabajo se han desarrollado y puesto a punto un conjunto de programas de ordenador escritos en FORTRAN 77, capaces de realizar todo el proceso de cálculo inherente al trabajo desarrollado. Los programas están optimizados para entornos multiprocesadores, tanto bajo paradigmas data-parallel como message-passing. Como librería de distribución de datos se ha empleado el Power Fortran Accelerator de SGI, y como paso de mensajes se ha empleado la librería estándar MPI. Todos los códigos serán disponibles en breve (gratuitamente) en la dirección internet zar.unizar.es.

El trabajo acometido es considerado por el autor como amplio y laborioso, pues la elaboración y desarrollo de programas requiere una cantidad de tiempo grande, especialmente si se persigue elaborar códigos precisos y eficientes. Tal y como dice uno de mis directores: “Saber calcular bien es una tarea difícil, conseguirlo realmente es mucho más difícil”.
Quedan pendientes numerosos trabajos que podrían acometerse en una segunda fase, como por ejemplo: introducción de anisotropía y visco-elasticidad en los modelos; introducción de librerías estandarizadas de elementos finitos o técnicas de resolución paralelas, que permitan dotar de una mayor versatilidad al paquete informático creado; adecuación del modelado de fuentes para la simulación de ondas superficiales…