TESIS Aplicaciones del Cómputo en Paralelo a la Modelación de Sistemas Terrestres”

UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO INSTITUTO DE GEOFÍSICA “Aplicaciones del Cómputo en Paralelo a la Modelación de Sistemas Terrestres” T E S I S QUE PARA OBTENER EL GRADO DE MAESTRO EN CIENCIAS DE LA TIERRA (MODELACIÓN DE SISTEMAS TERRESTRES) P R E S E N T A: MAT. ANTONIO CARRILLO LEDESMA DIRECTOR DE TESIS: DR. ISMAEL HERRERA REVILLA 2006 Dedico el presente trabajo con todo cariño a: Mi madre Alfonsina por haber hecho posible esto, después de pensar que yo era incorregible. Mi esposa Josefina por todo su apoyo y tiempo cedido para realizar la presente. Mi hijo José Antonio por mostrarme lo que realmente es la vida y ceder tanto de su tiempo de juegos conmigo para poder materializar este proyecto. Toda mi familia y amigos presentes y a los ya ausentes por mostrarme el camino en vida y así poder lograr este trabajo. Agradecimientos: Quiero agradecer a todas las personas que hicieron posible que el presente trabajo llegará a su fin: Al director de tesis: Ismael Herrera Revilla por todo el apoyo, tiempo y sus invaluables enseñanzas A los sinodales: Robert Yates Martín Díaz Viera Fabián García Nocetti Arón Jazcilevich Diamant por sus enseñanzas y comentarios A los investigadores: Alejandra Arciniega Ceballos Jorge Luís Ortega Arjona por todo el apoyo proporcionado Índice 1. Introducción 3 1.1. Antecedentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2. Métodos de Descomposición de Dominio . . . . . . . . . . . . . . 5 1.3. Objetivos de la Tesis . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3.1. Objetivos Generales . . . . . . . . . . . . . . . . . . . . . 8 1.3.2. Objetivos Particulares . . . . . . . . . . . . . . . . . . . . 8 1.4. Infraestructura Usada . . . . . . . . . . . . . . . . . . . . . . . . 9 2. Sistemas Continuos y sus Modelos 11 2.1. LosModelos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.1. FísicaMicroscópica y FísicaMacroscópica . . . . . . . . . 11 2.2. Cinemática de losModelos de Sistemas Continuos . . . . . . . . 12 2.2.1. Propiedades Intensivas y sus Representaciones . . . . . . 14 2.2.2. Propiedades Extensivas . . . . . . . . . . . . . . . . . . . 16 2.2.3. Balance de Propiedades Extensivas e Intensivas . . . . . . 17 2.3. Ejemplos deModelos . . . . . . . . . . . . . . . . . . . . . . . . . 20 3. Ecuaciones Diferenciales Parciales 23 3.1. Clasificación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.1. Condiciones Iniciales y de Frontera . . . . . . . . . . . . . 26 3.1.2. Modelos Completos . . . . . . . . . . . . . . . . . . . . . 27 3.2. Análisis Funcional y Problemas Variacionales . . . . . . . . . . . 29 3.2.1. Espacios de Sobolev . . . . . . . . . . . . . . . . . . . . . 29 3.2.2. Formulas de Green y Problemas Adjuntos . . . . . . . . . 33 3.2.3. Problemas Variacionales con Valor en la Frontera . . . . . 37 4. El Método Galerkin y el Método de Elemento Finito 41 4.1. Método Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2. Método de Elemento Finito . . . . . . . . . . . . . . . . . . . . . 44 4.2.1. Discretización Usando Rectángulos . . . . . . . . . . . . . 47 5. Solución de Grandes Sistemas de Ecuaciones 52 5.1. Métodos Directos . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.2. Métodos Iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.3. Precondicionadores . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.3.1. Gradiente Conjugado Precondicionado . . . . . . . . . . . 61 5.3.2. Precondicionador a Posteriori . . . . . . . . . . . . . . . . 63 5.3.3. Precondicionador a Priori . . . . . . . . . . . . . . . . . . 66 6. Métodos de Descomposición de Dominio (DDM) 69 6.1. Método de Schwarz . . . . . . . . . . . . . . . . . . . . . . . . . . 70 6.2. Método de Subestructuración . . . . . . . . . . . . . . . . . . . . 74 6.2.1. Precondicionador Derivado de laMatriz de Rigidez . . . . 82 1 7. El Cómputo en Paralelo 87 7.1. Arquitecturas de Software y Hardware . . . . . . . . . . . . . . . 87 7.1.1. Clasificación de Flynn . . . . . . . . . . . . . . . . . . . . 87 7.1.2. Categorías de Computadoras Paralelas . . . . . . . . . . . 89 7.2. Métricas de Desempeño . . . . . . . . . . . . . . . . . . . . . . . 94 7.3. Cómputo Paralelo para Sistemas Continuos . . . . . . . . . . . . 96 8. Implementación Computacional Secuencial y Paralela de DDM103 8.1. El Operador de Laplace y la Ecuación de Poisson . . . . . . . . . 104 8.2. Método del Elemento Finito Secuencial . . . . . . . . . . . . . . . 106 8.3. Método de Subestructuración Secuencial . . . . . . . . . . . . . . 108 8.4. Método de Subestructuración en Paralelo . . . . . . . . . . . . . 112 8.5. Método de Subestructuración en Paralelo Precondicionado . . . . 116 9. Análisis de Rendimiento y Conclusiones 119 9.1. Análisis de Comunicaciones . . . . . . . . . . . . . . . . . . . . . 119 9.2. Afectación del Rendimiento al Aumentar el Número de Subdominios en la Descomposición. . . . . . . . . . . . . . . . . . . . . 120 9.3. Descomposición Óptima para un Equipo Paralelo Dado. . . . . . 121 9.4. Consideraciones para Aumentar el Rendimiento . . . . . . . . . . 123 9.5. Conclusiones Generales . . . . . . . . . . . . . . . . . . . . . . . . 125 9.6. Trabajo Futuro . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 10.Apéndice 127 10.1.Nociones de Algebra Lineal . . . . . . . . . . . . . . . . . . . . . 127 10.2. σ-Algebra y EspaciosMedibles . . . . . . . . . . . . . . . . . . . 128 10.3. Espacios Lp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 10.4.Distribuciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 11.Bibliografía 134 2 1. Introducción 1.1. Antecedentes La necesidad de entender su entorno y anticiparse a los acontecimientos tiene raíces muy profundas en el ser humano. Desde la prehistoria, el hombre trató de predecir a la naturaleza, pues de ella dependía su supervivencia, para lo cual inicialmente nuestros antepasados utilizaron a la brujería, así como el pensamiento mágico y el religioso. Sin embargo, el medio más efectivo para predecir el comportamiento de la naturaleza es el método científico (o su antecesor el método empírico) y es por eso que este anhelo humano ancestral, a través de la historia, ha sido motor de la ciencia. La maduración y el progreso de la predicción científica es, sin duda, el resultado del avance general de la ciencia, pero además ha habido elementos catalizadores esenciales sin los cuales esto no hubiera sido posible. La predicción científica, además de ser científica, es matemática y computacional. En la actualidad, cuando deseamos predecir el comportamiento de un sistema, los conocimientos científicos y tecnológicos se integran en modelos matemáticos los cuales se convierten en programas de cómputo que son ejecutados por las computadoras tanto secuenciales como paralelas (entenderemos por una arquitectura paralela a un conjunto de procesadores interconectados capaces de cooperar en la solución de un problema). Aunque el sólo hecho de poseer la capacidad de predicción científica nos llena de satisfacción y orgullo, sin embargo, aún más trascendente es el hecho de que ella también es la base de una gran parte del extraordinario progreso material que la humanidad ha experimentado en épocas recientes. En efecto, nuestra facultad para predecir es una herramienta muy poderosa de la ingeniería, de la tecnología y de la ciencia misma, la cual, entre otras muchas cosas, nos ha permitido ampliar la disponibilidad de los recursos naturales y utilizarlos con mayor eficiencia. Para resolver, por ejemplo, algún problema del área de ciencias e ingeniería (el movimiento de un fluido libre o de un fluido en un medio poroso) lo primero que debemos de hacer es formular un modelo matemático del problema a tratar. Esto se logra a través de los fundamentos de la física macroscópica los cuales son proporcionados por la ‘teoría de los medios continuos’ [13]. Con base en ella se introduce una formulación clara, general y sencilla de los modelos matemáticos de los sistemas continuos. La formulación es tan sencilla y tan general, que los modelos básicos de sistemas tan complicados y diversos como la atmósfera, los océanos, los yacimientos petroleros, o los geotérmicos, se derivan por medio de la aplicación repetida de una sola ecuación diferencial: ‘la ecuación diferencial de balance’ [8]. Dicha formulación también es muy clara, pues en el modelo general no hay ninguna ambigüedad; en particular, todas las variables y parámetros que intervienen en él, están definidos de manera unívoca. En realidad, este modelo general de los sistemas continuos constituye una realización extraordinaria de los paradigmas del pensamiento matemático. El descubrimiento del hecho de que 3 los modelos matemáticos de los sistemas continuos, independientemente de su naturaleza y propiedades intrínsecas, pueden formularse por medio de balances, cuya idea básica no difiere mucho de los balances de la contabilidad financiera, fue el resultado de un largo proceso de perfeccionamiento en el que concurrieron una multitud de mentes brillantes. Los modelos matemáticos de los sistemas continuos son ecuaciones diferenciales, las cuales son parciales (con valores iniciales y condiciones de frontera) para casi todos los sistemas de mayor interés en la ciencia y la ingeniería, o sistemas de tales ecuaciones. Salvo para los problemas más sencillos, no es posible obtener por métodos analíticos las soluciones de tales ecuaciones, que son las que permiten predecir el comportamiento de los sistemas. La capacidad para formular los modelos matemáticos de sistemas complicados y de gran diversidad, es sin duda una contribución fundamental para el avance de la ciencia y sus aplicaciones, tal contribución quedaría incompleta y, debido a ello, sería poco fecunda, si no se hubiera desarrollado simultáneamente su complemento esencial: los métodos matemáticos y la computación electrónica. En cambio, la diversidad y complejidad de problemas que pueden ser tratados con métodos numéricos y computacionales es impresionante. Los modelos de los sistemas continuos -es decir, sistemas físicos macroscópicos tales como los yacimientos petroleros, la atmósfera, los campos electromagnéticos, los océanos, los metalúrgicos, el aparato circulatorio de los seres humanos, la corteza terrestre, lo suelos y las cimentaciones, muchos sistemas ambientales, y muchos otros cuya enumeración ocuparía un espacio enormecontienen un gran número de grados libertad. Por ello, la solución numérica por los esquemas tradicionales tipo diferencias finitas y elemento finito generan una discretización del problema, la cual es usada para generar sistemas de ecuaciones algebraicos [8]. Estos sistemas algebraicos en general son de gran tamaño para problemas reales, al ser estos algoritmos secuenciales su implantación suele hacerse en equipos secuenciales y por ello no es posible resolver muchos problemas que involucren el uso de una gran cantidad de memoria. Actualmente para tratar de subsanar la limitante de procesar sólo en equipos secuenciales, se usan equipos paralelos para soportar algoritmos secuenciales mediante directivas de compilación, haciendo ineficiente su implantación en dichos equipos. La computación en paralelo es una técnica que nos permite distribuir una gran carga computacional entre muchos procesadores. Y es bien sabido que una de las mayores dificultades del procesamiento en paralelo es la coordinación de las actividades de los diferentes procesadores y el intercambio de información entre los mismos [21]. Al usar métodos de descomposición de dominio conjuntamente con el cómputo en paralelo (supercomputadoras, clusters o grids) nos permite atacar una gran variedad de problemas que sin estas técnicas sería imposible hacerlo de manera eficiente. 4 1.2. Métodos de Descomposición de Dominio Los métodos de descomposición de dominio introducen desde la formulación matemática del problema una separación natural de las tareas a realizar por el método y simplifican considerablemente la transmisión de información entre los subdominios [2]. En ellos, los sistemas físicos representados por su modelo son descompuestos en varios subdominios. Los métodos de descomposición de dominio permiten tratar los problemas de tamaño considerable, empleando algoritmos paralelos en computadoras secuenciales y/o paralelas. Esto es posible ya que cualquier método de descomposición de dominio se basa en la suposición de que dado un dominio computacional Ω, este se puede particionar en subdominios Ωi,i = 1, 2, ..., E entre los cuales puede o no existir traslape [1] y [2]. Entonces, el problema es reformulado en términos de cada subdominio (empleando algún método del discretización) obteniendo una familia de subproblemas de tamaño reducido independientes en principio entre si, y que están acoplados a través de la solución en la interfaz de los subdominios que es desconocida. De esta manera, podemos clasificar de forma burda a los métodos de descomposición de dominio, como aquellos en que: existe traslape entre los subdominios y en los que no existe traslape. A la primera clase pertenece el método de Schwarz (en el cual el tamaño del traslape es importante en la convergencia del método) y a los de la segunda clase pertenecen los métodos del tipo subestructuración (en el cual los subdominios sólo tienen en común los nodos de la frontera interior). Estos métodos son un paradigma natural usado por la comunidad de modeladores. Los sistemas físicos son descompuestos en dos o más subdominios contiguos basados en consideraciones fenomenológicas. Estas descomposiciones basadas en dominios físicos son reflejadas en la ingeniería de software del código correspondiente. Así, mediante los métodos de descomposición de dominio, la programación orientada a objetos (que nos permite dividir en niveles la semántica de los sistemas complejos tratando así con las partes, más manejables que el todo, permitiendo su extensión y un mantenimiento más sencillo) y esquemas de paralelización que usan el paso de mensajes, es posible construir aplicaciones que coadyuven a la solución de problemas concomitantes en ciencia e ingeniería. Esta metodología permite utilizar todas las capacidades del cómputo en paralelo (grids de decenas de clusters, cada uno con miles de procesadores interconectados por red con un creciente poder de cómputo medible en Peta Flops), así como el uso de una amplia memoria (ya sea distribuida y/o compartida del orden de Tera Bytes), permitiendo atacar una gran variedad de problemas que sin estas técnicas es imposible hacerlo de manera flexible y eficiente. Pero hay que notar que existe una amplia gama de problemas que nos interesan resolver que superan las capacidades de cómputo actuales, ya sea por el tiempo requerido para su solución, por el consumo excesivo de memoria o ambos. La lista de los métodos de descomposición de dominio y el tipo de problemas 5 que pueden ser atacados por estos, es grande y está en constante evolución [2], ya que se trata de encontrar un equilibrio entre la complejidad del método (aunada a la propia complejidad del modelo), la eficiencia en el consumo de los recursos computacionales y la precisión esperada en la solución encontrada por los diversos métodos y las arquitecturas paralelas en la que se implante. 1.3. Objetivos de la Tesis Uno de los grandes retos del área de cómputo científico es poder analizar a priori una serie de consideraciones dictadas por factores externos al problema de interés que repercuten directamente en la forma de solucionar el problema, estas consideraciones influirán de manera decisiva en la implementación computacional de la solución numérica. Algunas de estas consideraciones son: • Número de Procesadores Disponibles • Tamaño y Tipo de Partición del Dominio • Tiempo de Ejecución Predeterminado Siendo común que ellas interactúan entre si, de forma tal que normalmente el encargado de la implementación computacional de la solución numérica tiene además de las complicaciones técnicas propias de la solución, el conciliarlas con dichas consideraciones. Esto deja al implementador de la solución numérica con pocos grados de libertad para hacer de la aplicación computacional una herramienta eficiente y flexible que cumpla con los lineamientos establecidos a priori y permita también que esta sea adaptable a futuros cambios de especificaciones -algo común en ciencia e ingeniería-. Para tratar de tener una idea clara de cómo afectan a la aplicación computacional estos factores, a continuación detallamos algunas consideraciones acerca de cada uno de ellos. Número de Procesadores Disponibles El número de procesadores disponibles es una barrera en constante evolución pero las principales limitantes son económicas y tecnológicas, es común ahora contar con miles de procesadores interactuando de manera conjunta en grids, pero todo ese poder de cómputo creciente es aun insuficiente para las necesidades del cómputo científico. Por otro lado la gran mayoría de los proyectos científicos y tecnológicos no cuentan con recursos computacionales grandes, ya que este es un recurso costoso y por ende muy limitado y en la mayoría de los casos no es el adecuado a las necesidades propias del problema. Así, para poder atacar el problema de forma eficiente es necesario adaptarse al equipo de cómputo disponible y tratar de conocer de antemano los factores que mermaran el rendimiento de la implementación computacional para buscar alternativas que mejoren la eficiencia de la implementación. 6 Tamaño y Tipo de Partición del Dominio Normalmente cuando se plantea un problema de sistemas continuos la determinación del dominio y tipo de malla a usar en la solución del mismo esta sujeta a restricciones fenomenológicas, por ello la malla deberá de adaptarse de la mejor forma posible para tratar de capturar los rasgos esenciales del fenómeno estudiado. Pero a la hora de la implementación computacional es común hacer adecuaciones a esta, para que pueda ser soportada por el equipo de cómputo y su ejecución sea en un tiempo razonable. Esto puede ser un problema sobre todos al ser implementada la solución usando cómputo paralelo, ya que comúnmente una elección de una malla no homogénea ocasionará problemas de mal balanceo de carga de trabajo entre los procesadores utilizados en la implementación computacional. Tiempo de Ejecución Predeterminado Otro factor determinante es el tiempo de solución esperado de la implementación computacional del problema, esto es crítico en problemas de control en tiempo real, comunes en la ciencia e ingeniería actual. Por ello en estos casos es permisible el aumento en el número y capacidades del equipo de cómputo necesario en la implementación con el fin de lograr un tiempo de ejecución por debajo del máximo permisible dictado por las especificaciones propias del problema. Todos estos factores influirán en la versión final implementada en la solución computacional de un problema particular, debiendo ser todas ellas sopesadas antes de tomar una decisión en cuanto a las especificaciones que el programa de cómputo deberá satisfacer. Así, el objetivo de este trabajo es desarrollar una metodología en la cual se implemente de forma paralela el método numérico de descomposición de dominio de subestructuración precondicionado, explicando detalladamente los fundamentos matemáticos de la modelación de fenómenos de sistemas continuos por medio de ecuaciones diferenciales parciales, los fundamentos del método de descomposición de dominio al aplicarse a ecuaciones diferenciales parciales elípticas y las ventajas sobre otros métodos de solución numérica. Mostraremos como construir el modelo computacional, esto no sólo nos ayudará a demostrar que es factible la construcción del propio modelo computacional a partir del modelo matemático y numérico para la solución de problemas reales. Además, mostrará los alcances y limitaciones en el consumo de recursos computacionales y nos permitirá la evaluación de algunas de las variantes de los métodos numéricos con los que es posible implementar el modelo computacional y haremos algo de análisis de rendimiento sin llegar a ser exhaustivo esté. También se muestra el diseño, análisis y programación de la aplicación computacional del método de descomposición de dominio en forma serial como paralela, basada en el paradigma de programación orientado a objetos y como este paradigma ofrece una forma robusta y flexible para la representación de entidades abstractas y que permiten una modelación más versátil así como una forma más eficiente para la organización y reutilización de código. 7 1.3.1. Objetivos Generales Mostrar las bases de una metodología que se utiliza para aplicar el cómputo en paralelo a la modelación matemática y computacional de sistemas continuos, de forma flexible, escalable y eficiente. Mostrar los alcances y limitaciones de la metodología usando como herramientas de evaluación a los métodos de elemento finito secuencial, método de subestructuración secuencial y método de subestructuración paralelo. Mostrar los diversos esquemas de optimización de los recursos computacionales aplicables a un problema especifico. Para facilitar la comprensión de las ideas básicas, se ha tomado la ecuación de Poisson. Es un ejemplo sencillo, pero gobierna los modelos de muchos sistemas de la ingeniería y de la ciencia, entre ellos el flujo de agua subterránea a través de un acuífero isotrópico, homogéneo bajo condiciones de equilibrio y es muy usada en múltiples ramas de la física. Por ejemplo, gobierna la ecuación de la conducción de calor en un sólido bajo condiciones de equilibrio. 1.3.2. Objetivos Particulares Como objetivos particulares de este trabajo tenemos: Mostrar cómo aplicar la metodología para manejar problemas de gran tamaño (descomposición de malla fina). Mostrar cómo descomponer un dominio en un conjunto de subdominios que den una partición en la que el tiempo de cálculo sea mínimo para una configuración de hardware dada. Mostrar cuales son las posibles optimizaciones aplicables a una configuración de hardware dada. Mostrar que es posible trabajar problemas con una malla muy fina en un equipo paralelo pequeño. Para poder cumplir con los objetivos, primeramente empezaremos describiendo las bases de los sistemas continuos y sus modelos, para conocer algunas propiedades de las ecuaciones diferenciales parciales que son generadas en la elaboración del modelo matemático. Después veremos la forma de pasar del modelo matemático al modelo numérico (el cual no es muy cercano al modelo computacional) introduciendo los fundamentos del método Galerkin y la derivación a partir de este del método de elemento finito, para posteriormente describir dos métodos de descomposición de dominio (Schwarz y subestructuración). Mostraremos como construir el modelo computacional que dependerá fuertemente de la arquitectura de cómputo disponible y del lenguaje de programación 8 usado. Esto nos ayudará a demostrar que es factible la construcción del propio modelo computacional a partir del modelo matemático y numérico para la solución de problemas reales. Además conoceremos los alcances y limitaciones en el consumo de recursos computacionales y nos permitirá la evaluación de algunas de las variantes de los métodos numéricos con los que es posible implementar el modelo computacional y haremos algo de análisis de rendimiento sin llegar a ser exhaustivo este. Finalmente se exploran los alcances y limitaciones de cada uno de los métodos implementados (elemento finito secuencial, descomposición de dominio secuencial y paralelo) y se muestra como es posible optimizar los recursos computacionales con que se cuente. 1.4. Infraestructura Usada El modelo computacional generado, está contenido en un programa de cómputo bajo el paradigma de orientación a objetos, programado en el lenguaje C++ en su forma secuencial y en su forma paralela en C++ en conjunto con la interfaz de paso de mensajes (MPI) bajo el esquema de maestroesclavo. Una versión de estos programas está disponible en la página Web http://www.mmc.igeofcu.unam.mx/, bajo el rubro FEM y DDM. Para desarrollar estos códigos, se realizo una jerarquía de clases para cada uno de los distintos componentes del sistema de elemento finito como de descomposición de dominio (rehusando todo el código de elemento finito en este), permitiendo usarlos tanto en forma secuencial como paralela. Las pruebas de rendimiento de los distintos programas se realizaron en equipos secuenciales y paralelos (clusters) que están montados en el Instituto de Geofísica de la UNAM, en las pruebas de análisis de rendimiento se usaron para la parte secuencial el equipo: • Computadora Pentium IV HT a 2.8 GHtz con 1 GB de RAM corriendo bajo el sistema operativo Linux Debian Stable con el compilador g++ de GNU. Para la parte paralela se usaron los siguientes equipos: • Cluster homogéneo de 10 nodos duales Xeon a 2.8 GHtz con 1 GB de RAM por nodo, unidos mediante una red Ethernet de 1 Gb, corriendo bajo el sistema operativo Linux Debian Stable con el compilador mpiCC de MPI de GNU. A cargo de la Dra. Alejandra Arciniega Ceballos del departamento de Geomagnetismo y Exploración. • Cluster heterogéneo con el nodo maestro Pentium IV HT a 3.4 GHtz con 1 GB de RAM y 7 nodos esclavos Pentium IV HT a 2.8 GHtz con 0.5 GB de RAM por nodo, unidos mediante una red Ethernet de 100 Mb, corriendo bajo el sistema operativo Linux Debian Stable con el compilador mpiCC de MPI de GNU. A cargo del Dr. Ismael Herrera Revilla del departamento de Recursos Naturales. 9 También se realizaron algunas pruebas de rendimiento en otro cluster heterogéneo de 64 nodos (donde algunos nodos son duales y otros hasta con 4 procesadores) Intel Xeon de 64 bits interconectados mediante una red de baja latencia a 1 Gb, pero los resultados no son comparables con respecto a los clusters anteriores debido a las diferencias entre las arquitecturas. Hay que notar que, el paradigma de programación orientada a objetos sacrifica algo de eficiencia computacional por requerir mayor manejo de recursos computacionales al momento de la ejecución. Pero en contraste, permite mayor flexibilidad a la hora adaptar los códigos a nuevas especificaciones. Adicionalmente, disminuye notoriamente el tiempo invertido en el mantenimiento y búsqueda de errores dentro del código. Esto tiene especial interés cuando se piensa en la cantidad de meses invertidos en la programación comparado con los segundos consumidos en la ejecución del mismo. 10 2. Sistemas Continuos y sus Modelos Los fundamentos de la física macroscópica los proporciona la ‘teoría de los medios continuos’. En este capítulo, con base en ella se introduce una formulación clara, general y sencilla de los modelos matemáticos de los sistemas continuos. Esta formulación es tan sencilla y tan general, que los modelos básicos de sistemas tan complicados y diversos como la atmósfera, los océanos, los yacimientos petroleros, o los geotérmicos, se derivan por medio de la aplicación repetida de una sola ecuación diferencial: ‘la ecuación diferencial de balance’. Dicha formulación también es muy clara, pues en el modelo general no hay ninguna ambigüedad; en particular, todas las variables y parámetros que intervienen en él, están definidos de manera unívoca. En realidad, este modelo general de los sistemas continuos constituye una realización extraordinaria de los paradigmas del pensamiento matemático. El descubrimiento del hecho de que los modelos matemáticos de los sistemas continuos, independientemente de su naturaleza y propiedades intrínsecas, pueden formularse por medio de balances, cuya idea básica no difiere mucho de los balances de la contabilidad financiera, fue el resultado de un largo proceso de perfeccionamiento en el que concurrieron una multitud de mentes brillantes. 2.1. Los Modelos Un modelo de un sistema es un sustituto de cuyo comportamiento es posible derivar el correspondiente al sistema original. Los modelos matemáticos, en la actualidad, son los utilizados con mayor frecuencia y también los más versátiles. En las aplicaciones específicas están constituidos por programas de cómputo cuya aplicación y adaptación a cambios de las propiedades de los sistemas es relativamente fácil. También, sus bases y las metodologías que utilizan son de gran generalidad, por lo que es posible construirlos para situaciones y sistemas muy diversos. Los modelos matemáticos son entes en los que se integran los conocimientos científicos y tecnológicos, con los que se construyen programas de cómputo que se implementan con medios computacionales. En la actualidad, la simulación numérica permite estudiar sistemas complejos y fenómenos naturales que sería muy costoso, peligroso o incluso imposible de estudiar por experimentación directa. En esta perspectiva la significación de los modelos matemáticos en ciencias e ingeniería es clara, porqué la modelación matemática constituye el método más efectivo de predecir el comportamiento de los diversos sistemas de interés. En nuestro país, ellos son usados ampliamente en la industria petrolera, en las ciencias y la ingeniería del agua y en muchas otras. 2.1.1. Física Microscópica y Física Macroscópica La materia, cuando se le observa en el ámbito ultramicroscópico, está formada por moléculas y átomos. Estos a su vez, por partículas aún más pequeñas como los protones, neutrones y electrones. La predicción del comportamiento de 11 estas partículas es el objeto de estudio de la mecánica cuántica y la física nuclear. Sin embargo, cuando deseamos predecir el comportamiento de sistemas tan grandes como la atmósfera o un yacimiento petrolero, los cuales están formados por un número extraordinariamente grande de moléculas y átomos, su estudio resulta inaccesible con esos métodos y en cambio el enfoque macroscópico es apropiado. Por eso en lo que sigue distinguiremos dos enfoques para el estudio de la materia y su movimiento. El primero -el de las moléculas, los átomos y las partículas elementales- es el enfoque microscópico y el segundo es el enfoque macroscópico. Al estudio de la materia con el enfoque macroscópico, se le llama física macroscópica y sus bases teóricas las proporciona la mecánica de los medios continuos. Cuando se estudia la materia con este último enfoque, se considera que los cuerpos llenan el espacio que ocupan, es decir que no tienen huecos, que es la forma en que los vemos sin el auxilio de un microscopio. Por ejemplo, el agua llena todo el espacio del recipiente donde está contenida. Este enfoque macroscópico está presente en la física clásica. La ciencia ha avanzado y ahora sabemos que la materia está llena de huecos, que nuestros sentidos no perciben y que la energía también está cuantizada. A pesar de que estos dos enfoques para el análisis de los sistemas físicos, el microscópico y el macroscópico, parecen a primera vista conceptualmente contradictorios, ambos son compatibles, y complementarios, y es posible establecer la relación entre ellos utilizando a la mecánica estadística. 2.2. Cinemática de los Modelos de Sistemas Continuos En la teoría de los sistemas continuos, los cuerpos llenan todo el espacio que ocupan. Y en cada punto del espacio físico hay una y solamente una partícula. Así, definimos como sistema continuo a un conjunto de partículas. Aún más, dicho conjunto es un subconjunto del espacio Euclidiano tridimensional. Un cuerpo es un subconjunto de partículas que en cualquier instante dado ocupa un dominio, en el sentido matemático, del espacio físico; es decir, del espacio Euclidiano tridimensional. Denotaremos por B(t) a la región ocupada por el cuerpo B, en el tiempo t, donde t puede ser cualquier número real. Figura 1: Representación del movimiento de partículas de un cuerpo B, para un tiempo dado. 12 Frecuentemente, sin embargo, nuestro interés de estudio se limitará a un intervalo finito de tiempo. Dado un cuerpo B, todo subdominio ˜ B ⊂ B, constituye a su vez otro cuerpo; en tal caso, se dice que ˜ B ⊂ B es un subcuerpo de B. De acuerdo con lo mencionado antes, una hipótesis básica de la teoría de los sistemas continuos es que en cualquier tiempo t ∈ (−∞, ∞) y en cada punto x ∈ B de la región ocupada por el cuerpo, hay una y sólo una partícula del cuerpo. Como en nuestra revisión se incluye no solamente la estática (es decir, los cuerpos en reposo), sino también la dinámica (es decir, los cuerpos en movimiento), un primer problema de la cinemática de los sistemas continuos consiste en establecer un procedimiento para identificar a las partículas cuando están en movimiento en el espacio físico. Sea X ∈ B, una partícula y p(X, t) el vector de la posición que ocupa, en el espacio físico, dicha partícula en el instante t. Una forma, pero no la única, de identificar la partícula X es asociándole la posición que ocupa en un instante determinado. Tomaremos en particular el tiempo t = 0, en tal caso p(X, 0) ≡ X. A las coordenadas del vector X ≡ (x1, x2, x3), se les llama las coordenadas materiales de la partícula. En este caso, las coordenadas materiales de una partícula son las coordenadas del punto del espacio físico que ocupaba la partícula en el tiempo inicial, t = 0. Desde luego, el tiempo inicial puede ser cualquier otro, si así se desea. Sea B el dominio ocupado por un cuerpo en el tiempo inicial, entonces X ∈ B si y solamente si la partícula X es del cuerpo. Es decir, B caracteriza al cuerpo. Sin embargo, debido al movimiento, la región ocupada por el mismo cambia con el tiempo y será denotada por B(t). Formalmente, para cualquier t ∈ (−∞, ∞), B(t) se define por B(t) ≡ © x ∈ R3 | ∃X ∈ B tal que x = p(X, t) ª (1) el vector posición p(X, t) es función del vector tridimensional X y del tiempo. Si fijamos el tiempo t, p(X, t) define una transformación del espacio Euclidiano R3 en si mismo y la Ec. (1) es equivalente a B(t) = p(B, t). Una notación utilizada para representar esta familia de funciones es p(•, t). De acuerdo a la hipótesis de los sistemas continuos: En cualquier tiempo t ∈ (−∞, ∞) y en cada punto x ∈ B de la región ocupada por el cuerpo hay una y sólo una partícula del cuerpo B para cada t fijo. Es decir, p(•, t) es una función biunívoca, por lo que existe la función inversa p−1(•, t). Si se fija la partícula X en la función p(X, t) y se varía el tiempo t, se obtiene su trayectoria. Esto permite obtener la velocidad de cualquier partícula, la cual es un concepto central en la descripción del movimiento. Ella se define como la derivada con respecto al tiempo de la posición cuando la partícula se mantiene fija. Es decir, es la derivada parcial con respecto al tiempo de la función de posición p(X, t). Por lo mismo, la velocidad como función de las coordenadas materiales de las partículas, está dada por V¯ (X, t) ≡ ∂p ∂t (X, t). (2) 13 2.2.1. Propiedades Intensivas y sus Representaciones En lo que sigue consideraremos funciones definidas para cada tiempo, en cada una de las partículas de un sistema continuo. A tales funciones se les llama ‘propiedades intensivas’. Las propiedades intensivas pueden ser funciones escalares o funciones vectoriales. Por ejemplo, la velocidad, definida por la Ec. (2), es una función vectorial que depende de la partícula X y del tiempo t. Una propiedad intensiva con valores vectoriales es equivalente a tres escalares, correspondientes a cada una de sus tres componentes. Hay dos formas de representar a las propiedades intensivas: la representación Euleriana y la representación Lagrangiana. Los nombres son en honor a los matemáticos Leonard Euler (1707-1783) y Joseph Louis Lagrange (1736-1813), respectivamente. Frecuentemente, el punto de vista Lagrangiano es utilizado en el estudio de los sólidos, mientras que el Euleriano se usa más en el estudio de los fluidos. Considere una propiedad intensiva escalar, la cual en el tiempo t toma el valor φ(X, t) en la partícula X. Entonces, de esta manera se define una función φ : B → R1, para cada t ∈ (−∞, ∞) a la que se denomina representación Lagrangiana de la propiedad intensiva considerada. Ahora, sea ψ(x, t) el valor que toma esa propiedad en la partícula que ocupa la posición x, en el tiempo t. En este caso, para cada t ∈ (−∞, ∞) se define una función ψ : B(t) → R1 a la cual se denomina representación Euleriana de la función considerada. Estas dos representaciones de una misma propiedad están relacionadas por la siguiente identidad φ(X, t) ≡ ψ(p(X, t), t). (3) Nótese que, aunque ambas representaciones satisfacen la Ec. (3), las funciones φ(X, t) y ψ(x, t) no son idénticas. Sus argumentos X y x son vectores tridimensionales (es decir, puntos de R3); sin embargo, si tomamos X = x, en general φ(X, t) 6= ψ(X, t). (4) La expresión de la velocidad de una partícula dada por la Ec. (2), define a su representación Lagrangiana, por lo que utilizando la Ec. (3) es claro que ∂p ∂t (X, t) =V ¯ (X, t) ≡v ¯ (p(X, t), t) (5) donde v ¯ (x, t) es la representación Euleriana de la velocidad. Por lo mismo v¯ (x, t) ≡V ¯ (p−1(x, t), t). (6) Esta ecuación tiene la interpretación de que la velocidad en el punto x del espacio físico, es igual a la velocidad de la partícula que pasa por dicho punto en el instante t. La Ec. (6) es un caso particular de la relación ψ(x, t) ≡ φ(p−1(x, t), t) de validez general, la cual es otra forma de expresar la relación de la Ec. (3) que existe entre las dos representaciones de una misma propiedad intensiva. 14 La derivada parcial con respecto al tiempo de la representación Lagrangiana φ(X, t) de una propiedad intensiva, de acuerdo a la definición de la derivada parcial de una función, es la tasa de cambio con respecto al tiempo que ocurre en una partícula fija. Es decir, si nos montamos en una partícula y medimos a la propiedad intensiva y luego los valores así obtenidos los derivamos con respecto al tiempo, el resultado final es ∂φ(X,t) ∂t . En cambio, si ψ(x, t) es la representación Euleriana de esa misma propiedad, entonces ∂ψ(x,t) ∂t es simplemente la tasa de cambio con respecto al tiempo que ocurre en un punto fijo en el espacio. Tiene interés evaluar la tasa de cambio con respecto al tiempo que ocurre en una partícula fija, cuando se usa la representación Euleriana. 2.2.2. Propiedades Extensivas En la sección anterior se consideraron funciones definidas en las partículas de un cuerpo, más precisamente, funciones que hacen corresponder a cada partícula y cada tiempo un número real, o un vector del espacio Euclidiano tridimensional R3. En ésta, en cambio, empezaremos por considerar funciones que a cada cuerpo B de un sistema continuo, y a cada tiempo t le asocia un número real o un vector de R3. A una función de este tipo E(B, t) se le llama ‘propiedad extensiva’ cuando esta dada por una integral E(B, t) ≡ Z B(t) ψ(x, t)dx¯ . (14) Observe que, en tal caso, el integrando define una función ψ(x, t) y por lo mismo, una propiedad intensiva. En particular, la función ψ(x, t) es la representación Euleriana de esa propiedad intensiva. Además, la Ec. (14) establece una correspondencia biunívoca entre las propiedades extensivas y las intensivas, porqué dada la representación Eulereana ψ(x, t) de cualquier propiedad intensiva, su integral sobre el dominio ocupado por cualquier cuerpo, define una propiedad extensiva. Finalmente, la notación empleada en la Ec. (14) es muy explicita, pues ahí se ha escrito E(B, t) para enfatizar que el valor de la propiedad extensiva corresponde al cuerpo B. Sin embargo, en lo que sucesivo, se simplificara la notación omitiendo el símbolo B es decir, se escribirá E(t) en vez de E(B, t). Hay diferentes formas de definir a las propiedades intensivas. Como aquí lo hemos hecho, es por unidad de volumen. Sin embargo, es frecuente que se le defina por unidad de masa véase [15]. Es fácil ver que la propiedad intensiva por unidad de volumen es igual a la propiedad intensiva por unidad de masa multiplicada por la densidad de masa (es decir, masa por unidad de volumen), por lo que es fácil pasar de un concepto al otro, utilizando la densidad de masa. Sin embargo, una ventaja de utilizar a las propiedades intensivas por unidad de volumen, en lugar de las propiedades intensivas por unidad de masa, es que la correspondencia entre las propiedades extensivas y las intensivas es más directa: dada una propiedad extensiva, la propiedad intensiva que le corresponde es la función que aparece como integrando, cuando aquélla se expresa como una integral de volumen. Además, del cálculo se sabe que ψ(x, t) ≡ l´ım V ol→0 E(t) V ol = l´ım V ol→0 R B(t) ψ(ξ, t)dξ V ol . (15) La Ec. (15) proporciona un procedimiento efectivo para determinar las propiedades extensivas experimentalmente: se mide la propiedad extensiva en un volumen pequeño del sistema continuo de que se trate, se le divide entre le volumen y el cociente que se obtiene es una buena aproximación de la propiedad intensiva. El uso que haremos del concepto de propiedad extensiva es, desde luego, lógicamente consistente. En particular, cualquier propiedad que satisface las condiciones de la definición de propiedad extensiva establecidas antes es, por 16 ese hecho, una propiedad extensiva. Sin embargo, no todas las propiedades extensivas que se pueden obtener de esta manera son de interés en la mecánica de los medios continuos. Una razón básica por la que ellas son importantes es porqué el modelo general de los sistemas continuos se formula en términos de ecuaciones de balance de propiedades extensivas, como se verá más adelante. 2.2.3. Balance de Propiedades Extensivas e Intensivas Los modelos matemáticos de los sistemas continuos están constituidos por balances de propiedades extensivas. Por ejemplo, los modelos de transporte de solutos (los contaminantes transportados por corrientes superficiales o subterráneas, son un caso particular de estos procesos de transporte) se construyen haciendo el balance de la masa de soluto que hay en cualquier dominio del espacio físico. Aquí, el término balance se usa, esencialmente, en un sentido contable. En la contabilidad que se realiza para fines financieros o fiscales, la diferencia de las entradas menos las salidas nos da el aumento, o cambio, de capital. En forma similar, en la mecánica de los medios continuos se realiza, en cada cuerpo del sistema continuo, un balance de las propiedades extensivas en que se basa el modelo. Ecuación de Balance Global Para realizar tales balances es necesario, en primer lugar, identificar las causas por las que las propiedades extensivas pueden cambiar. Tomemos como ejemplo de propiedad extensiva a las existencias de maíz que hay en el país. La primera pregunta es: ¿qué causas pueden motivar su variación, o cambio, de esas existencias?. Un análisis sencillo nos muestra que dicha variación puede ser debida a que se produzca o se consuma. También a que se importe o se exporte por los límites del país (fronteras o litorales). Y con esto se agotan las causas posibles; es decir, esta lista es exhaustiva. Producción y consumo son términos similares, pero sus efectos tienen signos opuestos, que fácilmente se engloban en uno solo de esos conceptos. De hecho, si convenimos en que la producción puede ser negativa, entonces el consumo es una producción negativa. Una vez adoptada esta convención, ya no es necesario ocuparnos separadamente del consumo. En forma similar, la exportación es una importación negativa. Entonces, el incremento en las existencias ΔE en un período Δt queda dado por la ecuación ΔE = P + I (16) donde a la producción y a la importación, ambas con signo, se les ha representado por P y I respectivamente. Similarmente, en la mecánica de los medios continuos, la lista exhaustiva de las causas por las que una propiedad extensiva de cualquier cuerpo puede cambiar, contiene solamente dos motivos: i) Por producción en el interior del cuerpo; y ii) Por importación (es decir, transporte) a través de la frontera. 17 Esto conduce a la siguiente ecuación de “balance global”, de gran generalidad, para las propiedades extensivas dE dt (t) = Z B(t) g(x, t)dx + Z ∂B(t) q(x, t)dx + Z Σ(t) gΣ(x, t)dx. (17) Donde g(x, t) es la generación en el interior del cuerpo, con signo, de la propiedad extensiva correspondiente, por unidad de volumen, por unidad de tiempo. Además, en la Ec. (17) se ha tomado en cuenta la posibilidad de que haya producción concentrada en la superficie Σ(t), la cual está dada en esa ecuación por la última integral, donde gΣ(x, t) es la producción por unidad de área. Por otra parte q(x, t) es lo que se importa o transporta hacia el interior del cuerpo a través de la frontera del cuerpo ∂B(t), en otras palabras, es el flujo de la propiedad extensiva a través de la frontera del cuerpo, por unidad de área, por unidad de tiempo. Puede demostrarse, con base en hipótesis válidas en condiciones muy generales, que para cada tiempo t existe un campo vectorial τ (x, t) tal que q(x, t) ≡ τ (x, t) • n(x, t) (18) donde n(x, t) es normal exterior a ∂B(t). En vista de esta relación, la Ec. (17) de balance se puede escribir como dE dt (t) = Z B(t) g(x, t)dx + Z ∂B(t) τ (x, t) • n(x, t)dx + Z Σ(t) gΣ(x, t)dx. (19) La relación (19) se le conoce con el nombre de “ecuación general de balance global” y es la ecuación básica de los balances de los sistemas continuos. A la función g(x, t) se le denomina el generación interna y al campo vectorial τ (x, t) el campo de flujo. Condiciones de Balance Local Los modelos de los sistemas continuos están constituidos por las ecuaciones de balance correspondientes a una colección de propiedades extensivas. Así, a cada sistema continuo le corresponde una familia de propiedades extensivas, tal que, el modelo matemático del sistema está constituido por las condiciones de balance de cada una de las propiedades extensivas de dicha familia. Sin embargo, las propiedades extensivas mismas no se utilizan directamente en la formulación del modelo, en su lugar se usan las propiedades intensivas asociadas a cada una de ellas. Esto es posible porqué las ecuaciones de balance global son equivalentes a las llamadas condiciones de balance local, las cuales se expresan en términos de las propiedades intensivas correspondientes. Las condiciones de balance local son de dos clases: ’las ecuaciones diferenciales de balance local’ y ’las condiciones de salto’. Las primeras son ecuaciones diferenciales parciales, que se deben satisfacer en cada punto del espacio ocupado por el sistema continuo, y las segundas son ecuaciones algebraicas que las discontinuidades deben satisfacer donde ocurren; es decir, en cada punto de Σ. Cabe mencionar que las ecuaciones diferenciales 18 de balance local son de uso mucho más amplio que las condiciones de salto, pues estas últimas solamente se aplican cuando y donde hay discontinuidades, mientras que las primeras en todo punto del espacio ocupado por el sistema continuo. Una vez establecidas las ecuaciones diferenciales y de salto del balance local, e incorporada la información científica y tecnológica necesaria para completar el modelo (la cual por cierto se introduce a través de las llamadas ’ecuaciones constitutivas’), el problema matemático de desarrollar el modelo y derivar sus predicciones se transforma en uno correspondiente a la teoría de la ecuaciones diferenciales, generalmente parciales, y sus métodos numéricos. Las Ecuaciones de Balance Local En lo que sigue se supone que las propiedades intensivas pueden tener discontinuidades, de salto exclusivamente, a través de la superficie Σ(t). Se entiende por ‘discontinuidad de salto’, una en que el límite por ambos lados de Σ(t) existe, pero son diferentes. Se utilizará en lo que sigue los resultados matemáticos que se dan a continuación, ver [8]. Teorema 1 Para cada t > 0, sea B(t) ⊂ R3 el dominio ocupado por un cuerpo. Suponga que la ‘propiedad intensiva’ ψ(x, t) es de clase C1, excepto a través de la superficie Σ(t). Además, sean las funciones v(x, t) y vΣ(x, t) esta última definida para x∈ Σ(t) solamente, las velocidades de las partículas y la de Σ(t), respectivamente. Entonces d dt Z B(t) ψdx ≡ Z B(t) ½ ∂ψ ∂t + ∇ • (vψ) ¾ dx + Z Σ [(v − vΣ) ψ] • ndx. (20) Teorema 2 Considere un sistema continuo, entonces, la ’ecuación de balance global’ (19) se satisface para todo cuerpo del sistema continuo si y solamente si se cumplen las condiciones siguientes: i) La ecuación diferencial ∂ψ ∂t + ∇ • (vψ) = ∇ • τ + g (21) vale en todo punto x∈ R3, de la región ocupada por el sistema. ii) La ecuación [ψ (v − vΣ) − τ ] • n = gΣ (22) vale en todo punto x∈ Σ. A las ecuaciones (21) y (22), se les llama ’ecuación diferencial de balance local’ y ’condición de salto’, respectivamente. Desde luego, el caso más general que se estudiará se refiere a situaciones dinámicas; es decir, aquéllas en que las propiedades intensivas cambian con el tiempo. Sin embargo, los estados estacionarios de los sistemas continuos son de sumo interés. Por estado estacionario se entiende uno en que las propiedades intensivas son independientes del tiempo. En los estados estacionarios, además, 19 las superficies de discontinuidad Σ(t) se mantienen fijas (no se mueven). En este caso ∂ψ ∂t = 0 y vΣ = 0. Por lo mismo, para los estados estacionarios, la ecuación de balance local y la condición de salto se reducen a ∇ • (vψ) = ∇ • τ + g (23) que vale en todo punto x∈ R3 y [ψv − τ ] • n = gΣ (24) que se satisface en todo punto de la discontinuidad Σ(t) respectivamente. 2.3. Ejemplos de Modelos Una de las aplicaciones más sencillas de las condiciones de balance local es para formular restricciones en el movimiento. Aquí ilustramos este tipo de aplicaciones formulando condiciones que se deben cumplir localmente cuando un fluido es incompresible. La afirmación de que un fluido es incompresible significa que todo cuerpo conserva el volumen de fluido en su movimiento. Entonces, se consideraran dos casos: el de un ‘fluido libre’ y el de un ‘fluido en un medio poroso’. En el primer caso, el fluido llena completamente el espacio físico que ocupa el cuerpo, por lo que el volumen del fluido es igual al volumen del dominio que ocupa el cuerpo, así Vf (t) = Z B(t) dx (25) aquí, Vf (t) es el volumen del fluido y B(t) es el dominio del espacio físico (es decir, de R3) ocupado por el cuerpo. Observe que una forma más explicita de esta ecuación es Vf (t) = Z B(t) 1dx (26) porqué en la integral que aparece en la Ec. (25) el integrando es la función idénticamente 1. Comparando esta ecuación con la Ec. (14), vemos que el volumen del fluido es una propiedad extensiva y que la propiedad intensiva que le corresponde es ψ = 1. Además, la hipótesis de incompresibilidad implica dVf dt (t) = 0 (27) esta es el balance global de la Ec. (19), con g = gΣ = 0 y τ = 0, el cual a su vez es equivalente a las Ecs. (21) y (22). Tomando en cuenta además que ψ = 1, la Ec. (21) se reduce a ∇ • v = 0. (28) Esta es la bien conocida condición de incompresibilidad para un fluido libre Además, aplicando la Ec. (22) donde haya discontinuidades, se obtiene [v]•n = 0. Esto implica que si un fluido libre es incompresible, la velocidad de sus partículas es necesariamente continua. 20 El caso en que el fluido se encuentra en un ‘medio poroso’, es bastante diferente. Un medio poroso es un material sólido que tiene huecos distribuidos en toda su extensión, cuando los poros están llenos de un fluido, se dice que el medio poroso esta ‘saturado’. Esta situación es la de mayor interés en la práctica y es también la más estudiada. En muchos de los casos que ocurren en las aplicaciones el fluido es agua o petróleo. A la fracción del volumen del sistema, constituido por la ‘matriz sólida’ y los huecos, se le llama ‘porosidad’ y se le representara por φ, así φ(x, t) = l´ım V →0 Volumen de huecos Volumen total (29) aquí hemos escrito φ(x, t) para enfatizar que la porosidad generalmente es función tanto de la posición como del tiempo. Las variaciones con la posición pueden ser debidas, por ejemplo, a heterogeneidad del medio y los cambios con el tiempo a su elasticidad; es decir, los cambios de presión del fluido originan esfuerzos en los poros que los dilatan o los encogen. Cuando el medio esta saturado, el volumen del fluido Vf es igual al volumen de los huecos del dominio del espacio físico que ocupa, así Vf (t) = Z B(t) φ(x, t)dx. (30) En vista de esta ecuación, la propiedad intensiva asociada al volumen de fluido es la porosidad φ(x, t) por lo que la condición de incomprensibilidad del fluido contenido en un medio poroso, esta dada por la ecuación diferencial ∂φ ∂t + ∇ • (vφ) = 0. (31) Que la divergencia de la velocidad sea igual a cero en la Ec. (28) como condición para que un fluido en su movimiento libre conserve su volumen, es ampliamente conocida. Sin embargo, este no es el caso de la Ec. (31), como condición para la conservación del volumen de los cuerpos de fluido contenidos en un medio poroso. Finalmente, debe observarse que cualquier fluido incompresible satisface la Ec. (28) cuando se mueve en el espacio libre y la Ec. (31) cuando se mueve en un medio poroso. Cuando un fluido efectúa un movimiento en el que conserva su volumen, al movimiento se le llama ‘isocorico’. Es oportuno mencionar que si bien cierto que cuando un fluido tiene la propiedad de ser incompresible, todos sus movimientos son isocoricos, lo inverso no es cierto: un fluido compresible en ocasiones puede efectuar movimientos isocoricos. Por otra parte, cuando un fluido conserva su volumen en su movimiento satisface las condiciones de salto de Ec. (22), las cuales para este caso son [φ(v − vΣ)] • n = 0. (32) En aplicaciones a geohidrología y a ingeniería petrolera, las discontinuidades de la porosidad están asociadas a cambios en los estratos geológicos y por esta 21 razón están fijas en el espacio; así, vΣ = 0 y la Ec. (32) se reduce a [φv] • n = 0 (33) o, de otra manera φ+vn+ = φ−vn− . (34) Aquí, la componente normal de la velocidad es vn ≡ v • n y los subíndices más y menos se utilizan para denotar los límites por los lado más y menos de Σ, respectivamente. Al producto de la porosidad por la velocidad se le conoce con el nombre de velocidad de Darcy U , es decir U = φv (35) utilizándola, las Ecs. (33) y (34) obtenemos [U ] • n = 0 y U n+ = U n− (36) es decir, 1. La Ec. (34) es ampliamente utilizada en el estudio del agua subterránea (geohidrología). Ahí, es frecuente que la porosidad φ sea discontinua en la superficie de contacto entre dos estratos geológicos diferentes, pues generalmente los valores que toma esta propiedad dependen por lo que vn+ 6= vn− de cada estrato. En tal caso, φ+ 6= φ− necesariamente. Para más detalles de la forma y del desarrollo de algunos modelos usados en ciencias de la tierra, véase [8], [15], [23] y [13]. 22 3. Ecuaciones Diferenciales Parciales Cada una de las ecuaciones de balance da lugar a una ecuación diferencial parcial u ordinaria (en el caso en que el modelo depende de una sola variable independiente), la cual se complementa con las condiciones de salto, en el caso de los modelos discontinuos. Por lo mismo, los modelos de los sistemas continuos están constituidos por sistemas de ecuaciones diferenciales cuyo número es igual al número de propiedades intensivas que intervienen en la formulación del modelo básico. Los sistemas de ecuaciones diferenciales se clasifican en elípticas, hiperbólicas y parabólicas. Es necesario aclarar que esta clasificación no es exhaustiva; es decir, existen sistemas de ecuaciones diferenciales que no pertenecen a ninguna de estas categorías. Sin embargo, casi todos los modelos de sistemas continuos, en particular los que han recibido mayor atención hasta ahora, si están incluidos en alguna de estas categorías. 3.1. Clasificación Es importante clasificar a las ecuaciones diferenciales parciales y a los sistemas de tales ecuaciones, porqué muchas de sus propiedades son comunes a cada una de sus clases. Así, su clasificación es un instrumento para alcanzar el objetivo de unidad conceptual. La forma más general de abordar la clasificación de tales ecuaciones, es estudiando la clasificación de sistemas de ecuaciones. Sin embargo, aquí solamente abordaremos el caso de una ecuación diferencia de segundo orden, pero utilizando un método de análisis que es adecuado para extenderse a sistemas de ecuaciones. La forma general de un operador diferencial cuasi-lineal de segundo orden definido en Ω ⊂ R2 es Lu ≡ a(x, y) ∂2u ∂x2 + b(x, y) ∂2u ∂x∂y + c(x, y) ∂2u ∂y2 = F (x, y, u, ux, uy) (37) para una función u de variables independientes x e y. Nos restringiremos al caso en que a, b y c son funciones sólo de x e y y no funciones de u. Para la clasificación de las ecuaciones de segundo orden consideraremos una simplificación de la ecuación anterior en donde F (x, y, u, ux, uy) = 0 y los coeficientes a, b y c son funciones constantes, es decir a ∂2u ∂x2 + b ∂2u ∂x∂y + c ∂2u ∂y2 = 0 (38) en la cual, examinaremos los diferentes tipos de solución que se pueden obtener para diferentes elecciones de a, b y c. Entonces iniciando con una solución de la forma u(x, y) = f (mx + y) (39) para una función f de clase C2 y para una constante m, que deben ser determinadas según los requerimientos de la Ec. (38). Usando un apostrofe para 23 denotar la derivada de f con respecto de su argumento, las requeridas derivadas parciales de segundo orden de la Ec. (38) son ∂2u ∂x2 = m2f 00, ∂2u ∂x∂y = mf 00, ∂2u ∂y2 = f 00 (40) sustituyendo la ecuación anterior en la Ec. (38) obtenemos ¡ am2 + bm + c ¢ f 00 = 0 (41) de la cual podemos concluir que f 00 = 0 ó am2+bm+c = 0 ó ambas. En el caso de que f 00 = 0 obtenemos la solución f = f0 + mx + y, la cual es una función lineal de x e y y es expresada en términos de dos constantes arbitrarias, f0 y m. En el otro caso obtenemos am2 + bm + c = 0 (42) resolviendo esta ecuación cuadrática para m obtenemos las dos soluciones m1 = ¡ −b + √2 b2 − 4ac ¢ 2a , m2 = ¡ −b − √2 b2 − 4ac ¢ 2a (43) de donde es evidente la importancia de los coeficientes de la Ec. (38), ya que el signo del discriminante ¡ b2 − 4ac ¢ es crucial para determinar el número y tipo de soluciones de la Ec. (42). Así, tenemos tres casos a considerar: Caso I. ¡ b2 − 4ac ¢ > 0, Ecuación Hiperbólica. La Ec. (42) tiene dos soluciones reales distintas, m1 y m2. Así cualquier función de cualquiera de los dos argumentos m1x + y ó m2x + y resuelven a la Ec. (38). Por lo tanto la solución general de la Ec. (38) es u(x, y) = F (m1x + y) + G(m2x + y) (44) donde F y G son cualquier función de clase C2. Un ejemplo de este tipo de ecuaciones es la ecuación de onda, cuya ecuación canónica es ∂2u ∂x2 − ∂2u ∂t2 = 0. (45) Caso II. ¡ b2 − 4ac ¢ = 0, Ecuación Parabólica. Asumiendo que b 6= 0 y a 6= 0 (lo cual implica que c 6= 0). Entonces se tiene una sola raíz degenerada de la Ec. (42) con el valor de m1 = −b 2a que resuelve a la Ec. (38). Por lo tanto la solución general de la Ec. (38) es u(x, y) = F (m1x + y) + yG(m1x + y) (46) donde F y G son cualquier función de clase C2. Si b = 0 y a = 0, entonces la solución general es u(x, y) = F (x) + yG(x) (47) 24 la cual es análoga si b = 0 y c = 0. Un ejemplo de este tipo de ecuaciones es la ecuación de difusión o calor, cuya ecuación canónica es ∂2u ∂x2 − ∂u ∂t = 0. (48) Caso III. ¡ b2 − 4ac ¢ < 0, Ecuación Elíptica. La Ec. (42) tiene dos soluciones complejas m1 y m2 las cuales satisfacen que m2 es el conjugado complejo de m1, es decir, m2 = m∗1 . La solución general puede ser escrita en la forma u(x, y) = F (m1x + y) + G(m2x + y) (49) donde F y G son cualquier función de clase C2. Un ejemplo de este tipo de ecuaciones es la ecuación de Laplace, cuya ecuación canónica es ∂2u ∂x2 + ∂2u ∂y2 = 0. (50) Consideremos ahora el caso de un operador diferencial lineal de segundo orden definido en Ω ⊂ Rn y consideremos también la ecuación homogénea asociada a este operador Lu = 0 (52) además, sea x∈ Ω un punto del espacio Euclidiano y V (x) una vecindad de ese punto. Sea una función u definida en V (x) con la propiedad de que exista una variedad Σ de dimensión n −1 cerrada y orientada, tal que la función u satisface la Ec. (52) en V (x)\Σ. Se supone además que existe un vector unitario n que apunta en la dirección positiva (único) está definido en Σ. Además, la función u y sus derivadas de primer orden son continuas a través de Σ, mientras que los límites de las segundas derivadas de u existen por ambos lados de Σ. Sea x∈ Σ tal que ∙ 6= 0 (53) para alguna pareja i, j = 1, ..., n. Entonces decimos que la función u es una solución débil de esta ecuación en x. Teorema 3 Una condición necesaria para que existan soluciones débiles de la ecuación homogénea (52) en un punto x∈ Σ es que Xn i=1 Xn j=1 aij ninj = 0. (54) 25 Así, si definimos a la matriz A= (aij) y observamos que n • A • n = Xn i=1 Xn j=1 aij ninj (55) entonces podemos decir que: I) Cuando todos los eigenvalores de la matriz A son distintos de cero y además del mismo signo, entonces se dice que el operador es Elíptico. II) Cuando todos los eigenvalores de la matriz A son distintos de cero y además n − 1 de ellos tienen el mismo signo, entonces se dice que el operador es Hiperbólico. III) Cuando uno y sólo uno de los eigenvalores de la matriz A es igual a cero, entonces se dice que el operador es Parabólico. Para el caso en que n = 2, esta forma de clasificación coincide con la dada anteriormente. 3.1.1. Condiciones Iniciales y de Frontera Dado un problema concreto de ecuaciones en derivadas parciales sobre un dominio Ω, si la solución existe, esta no es única ya que generalmente este tiene un número infinito de soluciones. Para que el problema tenga una y sólo una solución es necesario imponer condiciones auxiliares apropiadas y estas son las condiciones iniciales y condiciones de frontera. En esta sección sólo se enuncian de manera general las condiciones iniciales y de frontera que son esenciales para definir un problema de ecuaciones diferenciales: A) Condiciones Iniciales Las condiciones iniciales expresan el valor de la función al tiempo inicial t = 0 (t puede ser fijada en cualquier valor) u(x, y, 0) = γ(x, y). (56) B) Condiciones de Frontera Las condiciones de frontera especifican los valores que la función u(x, y, t) o ∇u(x, y, t) tomarán en la frontera ∂Ω, siendo de tres tipos posibles: 1) Condiciones tipo Dirichlet Especifica los valores que la función u(x, y, t) toma en la frontera ∂Ω u(x, y, t) = γ(x, y). (57) 26 2) Condiciones tipo Neumann Aquí se conoce el valor de la derivada de la función u(x, y, t) con respecto a la normal n a lo largo de la frontera ∂Ω ∇u(x, y, t) • n = γ(x, y). (58) 3) Condiciones tipo Robin Está condición es una combinación de las dos anteriores α(x, y)u(x, y, t) + β(x, y)∇u(x, y, t) • n = g∂(x, y) (59) ∀ x, y ∈ ∂Ω. En un problema dado se debe prescribir las condiciones iniciales al problema y debe de existir alguno de los tipos de condiciones de frontera o combinación de ellas en ∂Ω. 3.1.2. Modelos Completos Los modelos de los sistemas continuos están constituidos por: • Una colección de propiedades intensivas o lo que es lo mismo, extensivas. • El conjunto de ecuaciones de balance local correspondientes (diferenciales y de salto). • Suficientes relaciones que liguen a las propiedades intensivas entre sí y que definan a g, τ y v en términos de estas, las cuales se conocen como leyes constitutivas. Una vez que se han planteado las ecuaciones que gobiernan al problema, las condiciones iniciales, de frontera y mencionado los procesos que intervienen de manera directa en el fenómeno estudiado, necesitamos que nuestro modelo sea completo. Decimos que el modelo de un sistema es completo si define un problema bien planteado. Un problema de valores iniciales y condiciones de frontera es bien planteado si cumple que: i) Existe una y sólo una solución y, ii) La solución depende de manera continua de las condiciones iniciales y de frontera del problema. Es decir, un modelo completo es aquél en el cual se incorporan condiciones iniciales y de frontera que definen conjuntamente con las ecuaciones diferenciales un problema bien planteado. 27 A las ecuaciones diferenciales definidas en Ω ⊂ Rn Δu = 0 (60) ∂2u ∂t2 − Δu = 0 ∂u ∂t − Δu = 0 se les conoce con los nombres de ecuación de Laplace, ecuación de onda y ecuación del calor, respectivamente. Cuando se considera la primera de estas ecuaciones, se entiende que u es una función del vector x ≡ (x1, ..., xn), mientras que cuando se considera cualquiera de las otras dos, u es una función del vector x ≡ (x1, ..., xn, t). Así, en estos últimos casos el número de variables independientes es n + 1 y los conceptos relativos a la clasificación y las demás nociones discutidas con anterioridad deben aplicarse haciendo la sustitución n → n+1 e identificando xn+1 = t. Ecuación de Laplace Para la ecuación de Laplace consideraremos condiciones del tipo Robin. En particular, condiciones de Dirichlet y condiciones de Neumann. Sin embargo, en este último caso, la solución no es única pues cualquier función constante satisface la ecuación de Laplace y también ∂u ∂n = g∂ con g∂ = 0. Ecuación de Onda Un problema general importante consiste en obtener la solución de la ecuación de onda, en el dominio del espacio-tiempo Ω×[0, t], que satisface para cada t ∈ (0, t] una condición de frontera de Robin en ∂Ω y las condiciones iniciales u(x, 0) = u0(x) y ∂u ∂t (x, 0) = v0(x), ∀x ∈ Ω (61) aquí u0(x) y v0(x) son dos funciones prescritas. El hecho de que para la ecuación de onda se prescriban los valores iniciales, de la función y su derivada con respecto al tiempo, es reminiscente de que en la mecánica de partículas se necesitan las posiciones y las velocidades iniciales para determinar el movimiento de un sistema de partículas. Ecuación de Calor También para la ecuación del calor un problema general importante consiste en obtener la solución de la ecuación de onda, en el dominio del espacio-tiempo Ω × [0, t], que satisface para cada t ∈ (0, t] una condición de frontera de Robin en y ciertas condiciones iniciales. Sin embargo, en este caso en ellas sólo se prescribe a la función u(x, 0) = u0(x), ∀x ∈ Ω. (62) 28 3.2. Análisis Funcional y Problemas Variacionales Restringiéndonos a problemas elípticos, una cuestión central en la teoría de problemas elípticos con valores en la frontera se relaciona con las condiciones bajo las cuales uno puede esperar que el problema tenga solución y esta es única, así como conocer la regularidad de la solución. Todas estas preguntas se pueden responder, pero se requiere utilizar herramientas de análisis funcional y de problemas variacionales con valor en la frontera. En esta sección tratamos de dar un panorama de esta teoría mostrando los resultados que la sustentan y las condiciones bajo las cuales un problema elíptico tiene solución y esta es única, además de conocer la regularidad de la solución, para mayor referencia de estos resultados ver [12], [18] y [3]. 3.2.1. Espacios de Sobolev En esta subsección detallaremos algunos resultados de los espacios de Sobolev sobre el conjunto de números reales, en estos espacios son sobre los cuales trabajaremos tanto para plantear el problema elíptico como para encontrar la solución al problema. Primeramente definiremos lo que entendemos por un espacio L2. Definición 4 Una función medible u(x) definida sobre Ω ⊂ Rn se dice que pertenece al espacio L2(Ω) si Z Ω |u(x)|2 dx < ∞ (63) es decir, es integrable. La definición de los espacios medibles, espacios Lp, distribuciones y derivadas de distribuciones están dados en el apéndice, estos resultados son la base para poder definir a los espacios de Sobolev. Para poder expresar de forma compacta derivadas parciales de orden m o menor, usaremos la definición siguiente. Definición 5 Sea Zn + el conjunto de todas las n-duplas de enteros no negativos, un miembro de Zn + se denota usualmente por α ó β (por ejemplo α = (α1, α2, ..., αn). Denotaremos por |α| la suma |α| = α1+α2+...+αn y por Dαu la derivada parcial Dαu = ∂|α|u ∂xα1 1 ∂xα2 2 ...∂xαn n (64) así, si |α| = m, entonces Dαu denota la m-ésima derivada parcial de u. Entonces haciendo uso de la definición anterior, definimos lo que entendemos por un espacio de Sobolev. Definición 6 El espacio de Sobolev de orden m, denotado por Hm(Ω), es definido como el espacio que consiste de todas las funciones en L2(Ω) que satisfacen Hm(Ω) = © u | Dαu ∈ L2(Ω) ∀α tal que |α| ≤ m ª . (65) 29 El producto escalar h•, •i de dos elementos u y v ∈ Hm(Ω) esta dado por hu, viHm = Z Ω X |α|≤m (Dαu) (Dαv) dx para u, v ∈ Hm (Ω) . (66) Nota: Es común que el espacio L2(Ω) sea denotado por H0(Ω). Un espacio completo con producto interior es llamado un espacio de Hilbert, un espacio normado y completo es llamado espacio de Banach. Y como todo producto interior define una norma, entonces todo espacio de Hilbert es un espacio de Banach. Definición 7 La norma k•kHm inducida a partir del producto interior h•, •iHm queda definida por kuk2 Hm = hu, uiHm = Z Ω X |α|≤m (Dαu)2 dx. (67) Ahora, con norma k•kHm, el espacio Hm(Ω) es un espacio de Hilbert, esto queda plasmado en el siguiente resultado. Teorema 8 El espacio Hm(Ω) con la norma k•kHm es un espacio de Hilbert. Ya que algunas de las propiedades de los espacios de Sobolev sólo son validas cuando la frontera del dominio es suficientemente suave. Para describir al conjunto donde los espacios de Sobolev están definidos, es común pedirle algunas propiedades y así definimos lo siguiente. Definición 9 Una función f definida sobre un conjunto Γ ⊂ Rn es llamada Lipschitz continua si existe una constante L > 0 tal que |f (x) − f (y)| ≤ L |x − y| ∀x, y ∈ Γ. (68) Notemos que una función Lipschitz continua es uniformemente continua. Definición 10 Entenderemos por un dominio al conjunto Ω ⊂ Rn que sea abierto y conexo. Sea Ω ⊂ Rn (n ≥ 2) un dominio con frontera ∂Ω, sea x0 ∈ ∂Ω y construyamos la bola abierta con centro en x0 y radio ε, i.e. B(x0, ε), entonces definiremos el sistema coordenado (ξ1, ..., ξn) tal que el segmento ∂Ω ∩ B(x0, ε) pueda expresarse como una función ξn = f (ξ1, ..., ξn−1) (69) entonces definimos. Definición 11 La frontera ∂Ω del dominio Ω es llamada de Lipschitz si f definida como en la Ec. (69) es una función Lipschitz continua. 30 El siguiente teorema resume las propiedades más importantes de los espacios de Sobolev Hm (Ω) . Teorema 12 Sea Hm(Ω) el espacio de Sobolev de orden m y sea Ω ⊂ Rn un dominio acotado con frontera Lipschitz. Entonces i) Hr(Ω) ⊂ Hm(Ω) si r ≥ m ii) Hm(Ω) es un espacio de Hilbert con respecto a la norma k•kHm iii) Hm(Ω) es la cerradura con respecto a la norma k•kHm del espacio C∞(Ω). De la parte iii) del teorema anterior, se puede hacer una importante interpretación: Para toda u ∈ Hm(Ω) es siempre posible encontrar una función infinitamente diferenciable f, tal que este arbitrariamente cerca de u en el sentido que ku − fkHm < ε (70) para algún ε > 0 dado. Cuando m = 0, se deduce la propiedad H0(Ω) = L2(Ω) a partir del teorema anterior. Corolario 13 El espacio L2(Ω) es la cerradura, con respecto a la norma L2, del espacio C∞(Ω). Otra propiedad, se tiene al considerar a cualquier miembro de u ∈ Hm(Ω), este puede ser identificado con una función en Cm(Ω), después de que posiblemente sean cambiados algunos valores sobre un conjunto de medida cero, esto queda plasmado en los dos siguientes resultados. Teorema 14 Sean X y Y dos espacios de Banach, con X ⊂ Y. Sea i : X → Y tal que i (u) = u. Si el espacio X tiene definida la norma k•kX y el espacio Y tiene definida la norma k•kY , decimos que X está inmersa continuamente en Y si ki (u)kY = kukY ≤ K kukX (71) para alguna constante K > 0. Teorema 15 (Inmersión de Sobolev) Sea Ω ⊂ Rn un dominio acotado con frontera ∂Ω de Lipschitz. Si (m − k) > n/2, entonces toda función en Hm(Ω) pertenece a Ck(Ω), es decir, hay un miembro que pertenece a Ck(Ω). Además, la inmersión Hm(Ω) ⊂ Ck(Ω) (72) es continua. 31 Traza de una Función en Hm (Ω) . Una parte fundamental en los problemas con valores en la frontera definidos sobre el dominio Ω, es definir de forma única los valores que tomará la función sobre la frontera ∂Ω, en este apartado veremos bajo que condiciones es posible tener definidos de forma única los valores en la frontera ∂Ω tal que podamos definir un operador tr (•) continuo que actué en Ω tal que tr (u) = u|∂Ω. El siguiente lema nos dice que el operador tr (•) es un operador lineal continuo de C1 ¡ Ω ¢ a C (∂Ω), con respecto a las normas k•kH1(Ω) y k•kL2(∂Ω) . Lema 16 Sea Ω un dominio con frontera ∂Ω de Lipschitz. La estimación ktr (u)kL2(∂Ω) ≤ C kukH1(Ω) (73) se satisface para toda función u ∈ C1 ¡ Ω ¢ , para alguna constante C > 0. Ahora, para el caso tr (•) : H1 (Ω) → L2 (∂Ω) , se tiene el siguiente teorema. Teorema 17 Sea Ω un dominio acotado en Rn con frontera ∂Ω de Lipschitz. Entonces: i) Existe un único operador lineal acotado tr (•) : H1 (Ω) → L2 (∂Ω) , tal que ktr (u)kL2(∂Ω) ≤ C kukH1(Ω) , (74) con la propiedad que si u ∈ C1 ¡ Ω ¢ , entonces tr (u) = u |∂Ω . ii) El rango de tr (•) es denso en L2 (∂Ω). El argumento anterior puede ser generalizado para los espacios Hm (Ω) , de hecho, cuando m > 1, entonces para toda u ∈ Hm (Ω) tenemos que Dαu ∈ H1 (Ω) para |α| ≤ m − 1, (75) por el teorema anterior, el valor de Dαu sobre la frontera está bien definido y pertenece a L2 (Ω) , es decir tr (Dαu) ∈ L2 (Ω) , |α| ≤ m − 1. (76) Además, si u es m-veces continuamente diferenciable, entonces Dαu es al menos continuamente diferenciable para |α| ≤ m − 1 y tr (Dαu) = (Dαu) |∂Ω . (77) Definición de los Espacios Hm 0 (Ω) . Los espacio Hm 0 (Ω) surgen comúnmente al trabajar con problemas con valor en la frontera y serán aquellos espacios que se nulifiquen en la frontera del dominio, es decir Definición 18 Definimos a los espacios Hm 0 (Ω) como la cerradura, en la norma de Sobolev k•kHm , del espacio Cm 0 (Ω) de funciones con derivadas continuas del orden menor que m, todas las cuales tienen soporte compacto en Ω, es decir Hm 0 (Ω) es formado al tomar la unión de Cm 0 (Ω) y de todos los límites de sucesiones de Cauchy en Cm 0 (Ω) que no pertenecen a Cm 0 (Ω) . 32 Las propiedades básicas de estos espacios están contenidas en el siguiente resultado. Teorema 19 Sea Ω un dominio acotado en Rn con frontera ∂Ω suficientemente suave y sea Hm 0 (Ω) la cerradura de C∞0 (Ω) en la norma k•kHm , entonces a) Hm 0 (Ω) es la cerradura de C∞0 (Ω) en la norma k•kHm ; b) Hm 0 (Ω) ⊂ Hm(Ω); c) Si u ∈ Hm(Ω) pertenece a Hm 0 (Ω), entonces Dαu = 0, sobre ∂Ω, |α| ≤ m − 1. (78) 3.2.2. Formulas de Green y Problemas Adjuntos Una cuestión central en la teoría de problemas elípticos con valores en la frontera se relaciona con las condiciones bajo las cuales uno puede esperar una única solución a problemas de la forma Lu = fΩ en Ω ⊂ Rn (79) Bou = g0 B1u = g1 ... Bm−1u = gm−1 ⎫⎪⎪⎪⎬ ⎪⎪⎪⎭ en ∂Ω donde L es un operador elíptico de orden 2m, de forma Lu = X |α|≤m (−1)|α| Dα ⎛ ⎝ X |β|≤m aαβ(x)Dβ u ⎞ ⎠, x ∈ Ω ⊂ Rn (80) donde los coeficientes aαβ son funciones de x suaves y satisfacen las condiciones para que la ecuación sea elíptica, el conjunto B0, B1, ..., Bm−1 de operadores de frontera son de la forma Bj u = X |α|≤qj b(j) α Dαu = gj (81) y constituyen un conjunto de condiciones de frontera que cubren a L. Los coeficientes b(j) α son asumidos como funciones suaves. En el caso de problemas de segundo orden la Ec. (81) puede expresarse como una sola condición de frontera Bu = Xn j=1 bj ∂u ∂xj + cu = g en ∂Ω. (82) Antes de poder ver las condiciones bajo las cuales se garantice la existencia y unicidad es necesario introducir el concepto de formula de Green asociada con el operador L∗, para ello definimos: 33 Definición 20 Con el operador dado como en la Ec. (80), denotaremos por L∗ al operador definido por L∗u = X |α|≤m (−1)|α| Dα ⎛ ⎝ X |β|≤m aβα(x)Dβ u ⎞ ⎠ (83) y nos referiremos a L∗ como el adjunto formal del operador L. La importancia del adjunto formal es que si aplicamos el teorema de Green (67) a la integral Z Ω vLudx (84) obtenemos Z Ω vLudx = Z Ω uL∗vdx + Z ∂Ω F (u, v)ds (85) en la cual F (u, v) representa términos de frontera que se nulifican al aplicar el teorema ya que la función v ∈ H1 0 (Ω). Si L = L∗; i.e. aαβ = aβα el operador es llamado de manera formal el auto-adjunto. En el caso de problemas de segundo orden, dos sucesivas aplicaciones del teorema de Green (67) y obtenemos, para i y j fijos − tal que L es formalmente el auto-ajunto si aji = aij . Para hacer el tratamiento más simple, restringiremos nuestra atención al problema homogéneo, es decir, en el cual g0, g1, ..., gm−1 = 0 (esta no es una restricción real, ya que se puede demostrar que cualquier problema no-homogéneo con condiciones de frontera puede convertirse en uno con condiciones de frontera homogéneo de una manera sistemática), asumiremos también que Ω es suave y la frontera ∂Ω de Ω es de clase C∞. 34 Así, en lo que resta de la sección, daremos los pasos necesarios para poder conocer bajo que condiciones el problema elíptico con valores en la frontera del tipo Lu = fΩ en Ω ⊂ Rn (89) Bou = 0 B1u = 0 ... Bm−1u = 0 ⎫⎪⎪⎪⎬ ⎪⎪⎪⎭ en ∂Ω donde el operador L y Bj estan dados como en (80) y (81), con s ≥ 2m tiene solución y esta es única. Para ello, necesitamos adoptar el lenguaje de la teoría de operadores lineales, algunos resultados clave de algebra lineal están detallados en el apéndice. Primeramente denotemos N (Bj) al espacio nulo del operador de frontera Bj : Hs (Ω) → L2 (Ω) , entonces N (Bj) = {u ∈ Hs (Ω) | Bj u = 0 en ∂Ω} (90) para j = 0, 1, 2, ..., m − 1. Adicionalmente definimos al dominio del operador L, como el espacio D(L) = Hs (Ω) ∩ N (B0) ∩ ... ∩ N (Bm−1) (91) = {u ∈ Hs (Ω) | Bj u = 0 en ∂Ω, j = 0, 1, .., m − 1} . Entonces el problema elíptico con valores en la frontera de la Ec. (89) con s ≥ 2m, puede reescribirse como, dado L : D(L) → Hs−2m (Ω) , hallar u que satisfaga Lu = fΩ en Ω. (92) Lo primero que hay que determinar es el conjunto de funciones fΩ en Hs−2m (Ω) para las cuales la ecuación anterior se satisface, i.e. debemos identificar el rango R(L) del operador L. Pero como nos interesa conocer bajo que condiciones la solución u es única, entonces podemos definir el núcleo N (L) del operador L como sigue N (L) = {u ∈ D(L) | Lu = 0} (93) = {u ∈ Hs (Ω) | Lu = 0 en Ω, Bj u = 0 en ∂Ω, j = 0, 1, .., m − 1} . Si el N (L) 6= {0} , entonces no hay una única solución, ya que si u0 es una solución, entonces u0 +w también es solución para cualquier w ∈ N (L), ya que L(u0 + w) = Lu0 + Lw = Lu0 = fΩ. (94) Así, los elementos del núcleo N (L) de L deberán ser excluidos del dominio D(L) del operador L, para poder asegurar la unicidad de la solución u. 35 Si ahora, introducimos el complemento ortogonal N (L)⊥ del núcleo N (L) del operador L con respecto al producto interior L2, definiéndolo como N (L)⊥ = {v ∈ D(L) | (v, w) = 0 ∀w ∈ N (L)} . (95) De esta forma tenemos que D(L) = N (L) ⊕ N (L)⊥ (96) i.e. para toda u ∈ D(L), u se escribe como u = v + w donde v ∈ N (L)⊥ y w ∈ N (L). Además N (L) ∩ N (L)⊥ = {0} . De forma similar, podemos definir los espacios anteriores para el problema adjunto L∗u = fΩ en Ω ⊂ Rn (97) B∗o u = 0 B∗1 u = 0 ... B∗m−1u = 0 ⎫⎪⎪⎪⎬ ⎪⎪⎪⎭ en ∂Ω y definimos D(L∗) = Hs (Ω) ∩ N (B∗0 ) ∩ ... ∩ N (B∗m−1) (98) = © u ∈ Hs (Ω) | B∗j u = 0 en ∂Ω, j = 0, 1, .., m − 1 ª . Entonces el problema elíptico con valores en la frontera de la Ec. (89) con s ≥ 2m, puede reescribirse como, dado L∗ : D(L∗) → Hs−2m (Ω) , hallar u que satisfaga L∗u = fΩ en Ω. (99) Definiendo para el operador L∗ N (L∗) = {u ∈ D(L∗) | L∗u = 0} (100) = © u ∈ Hs (Ω) | L∗u = 0 en Ω, B∗j u = 0 en ∂Ω, j = 0, 1, .., m − 1 ª . y N (L∗)⊥ = {v ∈ D(L∗) | (v, w)L2 = 0 ∀w ∈ N (L∗)} . (101) Así, con estas definiciones, es posible ver una cuestión fundamental, esta es, conocer bajo que condiciones el problema elíptico con valores en la frontera de la Ec. (89) con s ≥ 2m tiene solución y esta es única, esto queda resuelto en el siguiente teorema cuya demostración puede verse en [3] y [12]. Teorema 21 Considerando el problema elíptico con valores en la frontera de la Ec. (89) con s ≥ 2m definido sobre un dominio Ω acotado con frontera ∂Ω suave. Entonces 36 i) Existe al menos una solución si y sólo si f ∈ N (L∗)⊥, esto es, si (f, v)L2(Ω) = 0 ∀v ∈ N (L∗). (102) ii) Asumiendo que la solución u existe, esta es única si u ∈ N (L)⊥, esto es, si (u, w)L2(Ω) = 0 ∀w ∈ N (L). (103) iii) Si existe una única solución, entonces existe una única constante C > 0, independiente de u, tal que kukHs ≤ C kfkHs−2m . (104) Nota: i) El teorema afirma que el operador L es un operador suprayectivo de D(L) sobre el subespacio de funciones en Hs−2m que satisface (103). Además el operador L es inyectivo si el dominio es restringido al espacio de funciones que satisfagan a (102). ii) La parte (iii) del teorema puede interpretarse como un resultado de regularidad, en el sentido en que se muestra u ∈ Hs−2m (Ω) si f ∈ Hs (Ω) . (105) 3.2.3. Problemas Variacionales con Valor en la Frontera Restringiéndonos ahora en problemas elípticos de orden 2 (problemas de orden mayor pueden ser tratados de forma similar), reescribiremos este en su forma variacional. La formulación variacional es más débil que la formulación convencional ya que esta demanda menor suavidad de la solución u, sin embargo cualquier problema variacional con valores en la frontera corresponde a un problema con valor en la frontera y viceversa. Además, la formulación variacional facilita el tratamiento de los problemas al usar métodos numéricos de ecuaciones diferenciales parciales, en esta sección veremos algunos resultados clave como es la existencia y unicidad de la solución de este tipo de problemas, para mayores detalles, ver [3] y [12]. Si el operador L está definido por Lu = −∇ • a • ∇u + cu (106) con a una matriz positiva definida, simétrica y c ≥ 0, el problema queda escrito como −∇ • a • ∇u + cu = fΩ en Ω (107) u = g en ∂Ω. 37 Si multiplicamos a la ecuación −∇ • a • ∇u + cu = fΩ por v ∈ V = H1 0 (Ω), podemos reescribir el problema dado por la Ec. (79) de orden 2, haciendo uso de la forma bilineal a (•, •) y la funcional lineal l (•). Entonces entenderemos en el presente contexto un problema variacional con valores de frontera (VBVP) por uno de la forma: hallar una función u que pertenezca a un espacio de Hilbert V = H1 0 (Ω) y que satisfaga la ecuación a (u, v) = hf, vi (112) para toda función v ∈ V donde a (•, •) es una forma bilineal y l (•) es una funcional lineal. Definición 22 Sea V un espacio de Hilbert y sea k•kV la norma asociada a dicho espacio, decimos que una forma bilineal a (•, •) es continua si existe una constante M > 0 tal que |a (u, v)| ≤ M kukV kvkV ∀u, v ∈ V (113) y es V -elíptica si existe una constante α > 0 tal que a (v, v) ≥ α kvk2 V ∀v ∈ V (114) donde k•kV es la norma asociada al espacio V. Esto significa que una forma V − elíptica es una que siempre es no negativa y toma el valor de 0 sólo en el caso de que v = 0, i.e. es positiva definida. Notemos que el problema (107) definido en V = H1 0 (Ω) reescrito como el problema (112) genera una forma bilineal V -elíptica cuyo producto interior sobre V es simétrico y positivo definido ya que a (v, v) ≥ α kvk2 V > 0, ∀v ∈ V, v 6= 0 (115) 38 reescribiéndose el problema (112), en el cual debemos encontrar u ∈ V tal que a (u, v) = hf, vi − a (u0, v) (116) donde u0 = g en ∂Ω, para toda v ∈ V. Entonces, la cuestión fundamental, es conocer bajo que condiciones el problema anterior tiene solución y esta es única, el teorema de Lax-Milgram nos da las condiciones bajo las cuales el problema (107) reescrito como el problema (112) tiene solución y esta es única, esto queda plasmado en el siguiente resultado. Teorema 23 (Lax-Milgram) Sea V un espacio de Hilbert y sea a (•, •) : V × V → R una forma bilineal continua V -elíptica sobre V. Además, sea l (•) : V → R una funcional lineal continua sobre V. Entonces i) El VBVP de encontrar u ∈ V que satisfaga a (u, v) = hf, vi , ∀v ∈ V (117) tiene una y sólo una solución; ii) La solución depende continuamente de los datos, en el sentido de que kukV ≤ 1 α klkV ∗ (118) donde k•kV ∗ es la norma en el espacio dual V ∗ de V y α es la constante de la definición de V -elíptica. Más específicamente, considerando ahora V un subespacio cerrado de Hm(Ω) las condiciones para la existencia, unicidad y la dependencia continua de los datos queda de manifiesto en el siguiente resultado. Teorema 24 Sea V un subespacio cerrado de Hm(Ω), sea a (•, •) : V × V → R una forma bilineal continua V -elíptica sobre V y sea l (•) : V → R una funcional lineal continua sobre V. Sea P un subespacio cerrado de V tal que a (u + p, v + p) = a (u, v) ∀u, v ∈ V y p, p ∈ P. (119) También denotando por Q el subespacio de V consistente de las funciones ortogonales a P en la norma L2; tal que Q = ½ v ∈ V | Z Ω updx = 0 ∀p ∈ P ¾ , (120) y asumiendo que a (•, •) es Q-elíptica: existe una constante α > 0 tal que a (q, q) ≥ α kqk2 Q para q ∈ Q, (121) 39 la norma sobre Q será la misma que sobre V. Entonces i) Existe una única solución al problema de encontrar u ∈ Q tal que a (u, v) = hl, vi , ∀v ∈ V (122) si y sólo si las condiciones de compatibilidad hl, pi = 0 para p ∈ P (123) se satisfacen. ii) La solución u satisface kukQ ≤ α−1 klkQ∗ (124) (dependencia continua de los datos). Otro aspecto importante es la regularidad de la solución, si la solución u al VBVP de orden 2m con f ∈ Hs−2m(Ω) donde s ≥ 2m, entonces u pertenecerá a Hs(Ω) y esto queda de manifiesto en el siguiente resultado. Teorema 25 Sea Ω ⊂ Rn un dominio suave y sea u ∈ V la solución al VBVP a (u, v) = hf, vi, v ∈ V (125) donde V ⊂ Hm(Ω). Si f ∈ Hs−2m(Ω) con s ≥ 2m, entonces u ∈ Hs(Ω) y la estimación kukHs ≤ C kfkHs−2m (126) se satisface. 40 4. ElMétodo Galerkin y elMétodo de Elemento Finito En el presente capítulo se prestará atención a varios aspectos necesarios para encontrar la solución aproximada de problemas con valor en la frontera (VBVP). Ya que en general encontrar la solución a problemas con geometría diversa es difícil y en algunos casos imposible usando métodos analíticos. Por ello se explorará el método Galerkin para obtener la solución aproximada al problema de ecuaciones diferenciales parciales. El método de elemento finito está basado en la formulación del método Galerkin, tomando como punto de partida la formulación variacional. Pero hay que notar que existen una gama amplia de métodos para la solución de VBVP, algunos derivados del método Galerkin y otros que utilizan otro enfoque a los vistos en el presente trabajo, como los desarrollados en [16] conocidos como métodos de Trefftz-Herrera. En este capítulo se considera el VBVP de la forma Lu = fΩ en Ω (127) u = g en ∂Ω donde Lu = −∇ • a • ∇u + cu (128) con a una matriz positiva definida, simétrica y c ≥ 0, como un caso particular del operador elíptico definido por la Ec. (79) de orden 2, con Ω ⊂ R2 un dominio poligonal, es decir, Ω es un conjunto abierto acotado y conexo tal que ¯Ω es la unión de un número finito de polígonos. La sencillez del operador L nos permite facilitar la comprensión de muchas de las ideas básicas que se expondrán a continuación, pero tengamos en mente que esta es una ecuación que gobierna los modelos de muchos sistemas de la ciencia y la ingeniería, por ello es muy importante su solución. Y que al considerar al operador definido así, nos permite garantizar que la forma bilineal asociada será simétrica y definida positiva, lo cual como será visto posteriormente permite usar una serie de herramientas que sacan ventaja del operador, como es el uso de el método de gradiente conjugado usado para encontrar la solución del sistema algebraico de ecuaciones generado al buscar la solución aproximada al VBVP. 4.1. Método Galerkin Si multiplicamos a la ecuación −∇ • a • ∇u + cu = fΩ por v ∈ V = H1 0 (Ω), obtenemos −v ¡ ∇ • a • ∇u + cu ¢ = vfΩ (129) aplicando el teorema de Green (67) obtenemos la Ec. (86), que podemos reescribir como Z Ω ¡ ∇v • a • ∇u + cuv ¢ dx = Z Ω vfΩdx. (130) 41 Definiendo el operador bilineal a (u, v) = Z Ω ¡ ∇v • a • ∇u + cuv ¢ dx (131) y la funcional lineal l(v) = hf, vi = Z Ω vfΩdx (132) podemos reescribir el problema dado por la Ec. (127) de orden 2 en forma variacional, haciendo uso de la forma bilineal a (•, •) y la funcional lineal l (•). La idea básica detrás del método Galerkin es, considerando el VBVP, encontrar u ∈ V = H1 0 (Ω) que satisfaga a (u, v) = hf, vi ∀v ∈ V (133) donde V es un subespacio de un espacio de Hilbert H (por conveniencia nos restringiremos a espacios definidos sobre los números reales). El problema al tratar de resolver la Ec. (133) está en el hecho de que el espacio V es de dimensión infinita, por lo que resulta que en general no es posible encontrar el conjunto solución. En lugar de tener el problema en el espacio V, se supone que se tienen funciones linealmente independientes φ1, φ2, ..., φN en V y definimos el espacio V h a partir del subespacio dimensionalmente finito de V generado por las funciones φi, es decir, V h = Generado {φi}N i=1 , V h ⊂ V. (134) El índice h es un parámetro que estará entre 0 y 1, cuya magnitud da alguna indicación de cuan cerca V h esta de V, h se relaciona con la dimensión de V h. Y como el número N de las funciones base se escoge de manera que sea grande y haga que h sea pequeño. En el límite, cuando N → ∞, h → 0, se puede elegir {φi} tal que V h se aproxime a V de la forma que se detallará posteriormente. Después de definir el espacio V h, es posible trabajar con V h en lugar de V y encontrar una función uh que satisfaga a (uh, vh) = hf, vhi ∀vh ∈ V h. (135) Esta es la esencia del método Galerkin, notemos que uh y vh son sólo combinaciones lineales de las funciones base de V h, tales que uh = XN i=1 ciφi y vh = XN j=1 dj φj (136) donde vh es arbitraria, como los coeficientes de dj y sin perdida de generalidad podemos hacer vh = φj . Así, para encontrar la solución uh sustituimos las Ecs. (136) en la Ec. (135) y usando el hecho que a (•, •) es una forma bilineal y l (•) es una funcional lineal notemos que tanto Kij y Fj pueden ser evaluados, ya que φi, a (•, •) y l (•) son conocidas. Entonces el problema se reduce a resolver el sistema de ecuaciones lineales XN i=1 Kij ci − Fj, j = 1, 2, ...N (140) o más compactamente Ku = F (141) en la cual K y F son la matriz y el vector cuyas entradas son Kij y Fj . Una vez que el sistema es resuelto, la solución aproximada uh es encontrada. Notemos que el problema (127) definido en V h = H1 0 (Ω) reescrito como el problema (133) genera una forma bilineal V h-elíptica cuyo producto interior sobre V h es simétrico y positivo definido ya que a (vh, vh) ≥ α kvhk2 V h > 0, ∀vh ∈ V h, vh 6= 0 (142) reescribiéndose el problema (135) como el problema aproximado en el cual debemos encontrar uh ∈ V h ⊂ V tal que a (uh, vh) = hf, vhi − a (u0, vh) (143) donde u0 = g = 0 en ∂Ω, para toda vh ∈ V h, es decir Z Ω ¡ ∇vh • a • ∇uh + cuhvh ¢ dxdy = Z Ω fΩvhdxdy (144) para todo vh ∈ V h. Entonces, el problema (127) al aplicarle el método Galerkin obtenemos (130), el cual podemos reescribirlo como (144). Aplicando el teorema de Lax-Milgram (23) a este caso particular, tenemos que este tiene solución única y esta depende continuamente de los datos. El siguiente resultado nos da una condición suficiente para que la aproximación uh del método Galerkin converja a la solución u del problema dado por la Ec. (133), para más detalle véase [12] y [3]. Teorema 26 Sea V un subespacio cerrado de un espacio de Hilbert, y sea la forma bilineal a (•, •) : V h×V h → R continua V -elíptica y sea l (•) una funcional lineal acotada. Entonces existe una constante C, independiente de h, tal que ku − uhkV ≤ C inf vh∈V h ku − vhkV (145) 43 donde u es solución de (133) y uh es solución de (143), consecuentemente, una condición suficiente para que la aproximación uh del método Galerkin converge a la solución u del problema © dado por la Ec. (133) es que exista una familia V h ª de subespacios con la propiedad de que inf vh∈V h ku − vhkV → 0 cuando h → 0. (146) La derivación de algunos otros métodos tomando como base al método Galerkin pueden consultarse en [3], [15], [25] y [18]. 4.2. Método de Elemento Finito El método de elementos finitos provee una manera sistemática y simple de generar las funciones base en un dominio con geometría Ω poligonal. Lo que hace al método de elemento finito especialmente atractivo sobre otros métodos, es el hecho de que las funciones base son polinomios definidos por pedazos (elementos Ωi) que son no cero sólo en una pequeña parte de Ω, proporcionando a la vez una gran ventaja computacional al método ya que las matrices generadas resultan bandadas ahorrando memoria al implantarlas en una computadora. Así, partiendo del problema aproximado (144), se elegirá una familia de espacios V h(h ∈ (0, 1)) definido por el procedimiento de elementos finitos (descritos en las subsecciones siguientes en el caso de interpoladores lineales, para otros tipos de interpoladores, ver [15]), teniendo la propiedad de que V h se aproxima a V cuando h se aproxima a cero en un sentido apropiado, esto es, por supuesto una propiedad indispensable para la convergencia del método Galerkin. Mallado del dominio Para comenzar con el método, dividimos el dominio Ω ⊂ R2 en E subdominios o elementos Ωe llamados elementos finitos, tal que Ω = [E e=1 Ωe donde: • Cada Ωe es un polígono (rectángulo o triángulo) con interior no vacío (°Ωe 6= ∅). • Ωi ∩ Ωj = ∅ para cada i 6= j. • El diámetro Diam(Ωe) ≤ h para cada e = 1, 2, ..., E. • Los vértices de cada Ωe son llamados nodos, teniendo N de ellos. Funciones Base A continuación describiremos la manera de construir las funciones base usada por el método de elemento finito. En este procedimiento debemos tener en cuenta que las funciones base están definidas en un subespacio de V = H1 (Ω) para problemas de segundo orden que satisfacen las condiciones de frontera. Las funciones base deberán satisfacer las siguientes propiedades: 44 i) Las funciones base φi son acotadas y continuas, i.e φi ∈ C (Ωe) . ii) Existen 􀁣 funciones base por cada nodo del polígono Ωe, y cada función φi es no cero solo en los elementos contiguos conectados por el nodo i. iii) φi = 1 en cada i nodo del polígono Ωe y cero en los otros nodos. iv) La restricción φi a Ωe es un polinomio, i.e. φi ∈ Pk[Ωe] para alguna k ≥ 1 donde Pk[Ωe] es el espacio de polinomios de grado a lo más k sobre Ωe. Decimos que φi ∈ Pk[Ωe] es una base de funciones y por su construcción es evidente que estas pertenecen a H1 (Ω) . Al conjunto formado por todas las funciones base definidas para todo Ωe de Ω será el espacio Ph[k] de funciones base, i.e. Ph[k] = [E e=1 Pk[Ωe] estas formarán las funciones base globales. Solución aproximada Para encontrar la solución aproximada elegimos el espacio Ph[k] de funciones base, como el espacio de funciones lineales φi definidas por pedazos de grado menor o igual a k (en nuestro caso k = 1), entonces el espacio a trabajar es V h = Generado © φi ∈ Ph[k] | φi(x) = 0 en ∂Ω ª . (147) La solución aproximada de la Ec. (144) al problema dado por la Ec. (127) queda en términos fΩφj dxdy (150) entonces la matriz K ≡ [Kij ], los vectores u ≡ (u1, ..., uN ) y F ≡ (F1, ..., FN ) definen el sistema lineal (que es positivo definido) Ku = F (151) donde u será el vector solución a la Ec. (151) cuyos valores serán la solución al problema dado por la Ec. (144) que es la solución aproximada a la Ec. (127) en los nodos interiores de Ω. 45 Un Caso más General Sea el operador elíptico (caso simétrico) en el dominio Ω, y el operador definido por Lu = fΩ en Ω\Σ (152) u = g en ∂Ω [u]Σ = J0 [an • ∇u]Σ = J1 donde Lu = −∇ • a • ∇u + cu (153) conjuntamente con una partición Q = {Ω1, ..., ΩE} de Ω. Multiplicando por la función w obtenemos wLu = −w∇ • a • ∇u + cwu = wfΩ (154) entonces si w(x) es tal que [w] = 0 (es decir w es continua) y definimos a (u, w) = XE i=1 Z Ωi ¡ ∇u • a • ∇w + cwu ¢ dx (155) tal que a (u, w) define un producto interior sobre H1 (Ω) = H1 (Ω1) ⊕ H1 (Ω2) ⊕ .... ⊕ H1 (ΩE) . Sea u0(x) una función que satisface las condiciones de frontera y J0 una función que satisface las condiciones de salto, tal que i) u0(x) = g(x) en ∂Ω ii) [u0(x)]Σ = J0 y sea u(x) = u0(x) + v(x). Entonces u(x) satisface la Ec. (155) si y sólo si v(x) satisface a (u, w) = Z Ω wfΩdx − hu0, wi − Z Σ J1wds (157) para toda w tal que w(x) = 0 en ∂Ω. Sea {φi} una base de un subespacio de dimensión finita V h definido como V h = © φi | φi ∈ C1 (Ωi) , ∀i, φi = 0 en ∂Ω y φi ∈ C0 (Ω) ª . (158) 46 La solución por elementos finitos de (157) se obtiene al resolver el sistema lineal esta solución será la solución en los nodos interiores de Ω. 4.2.1. Discretización Usando Rectángulos Para resolver la Ec. (127), usando una discretización con rectángulos, primero dividimos el dominio Ω ⊂ R2 en Nx nodos horizontales por Ny nodos verticales, teniendo E = (Nx − 1)(Ny − 1) subdominios o elementos rectangulares Ωe tales que Ω = ∪E e=1Ωe y Ωi ∩ Ωj 6= ∅ si son adyacentes, con un total de N = NxNy nodos. Donde las funciones lineales definidas por pedazos en Ωe en nuestro caso serán polinomios de orden uno en cada variable separadamente y cuya restricción de φi a Ωe es φ(e) i la cual deberá ser ensamblada en la matriz de carga global que corresponda a la numeración de nodos locales del elemento Ωe con respecto a la numeración global de los elementos en Ω. De manera parecida, para cada Ωe de Ω se genera el vector de integrales (vector de carga local) Fj = Z Ωe fΩφ(e) j dxdy (166) con la estructura ⎡ ⎢⎢⎢⎣ F (e) 1 F (e) 2 F (e) 3 F (e) 4 ⎤ ⎥⎥⎥⎦ el cual también deberá ser ensamblado en el vector de carga global que corresponda a la numeración de nodos locales al elemento Ωe con respecto a la numeración global de los elementos de Ω. Montando los K(e) ij en la matriz K y los F (e) j en el vector F según la numeración de nodos global, se genera el sistema Kuh=F donde uh será el vector cuyos valores serán la solución aproximada a la Ec. (127) en los nodos interiores de Ω. La matriz K generada de esta forma, tiene una propiedad muy importante, es bandada y el ancho de banda es de 9 elementos, esto es muy útil al momento de soportar la matriz en memoria. teniendo en mente el simplificar los cálculos computacionales, se considera un elemento de referencia ˆΩ en los ejes coordenados (ε, η) cuyos vértices están el (−1, −1), (1, −1), (1, 1) y (−1, 1) respectivamente, en el cual mediante una 48 función afín será proyectado cualquier elemento rectangular Ωe cuyos vértices (x(e) 1 , y(e) 1 ), (x(e) 2 , y(e) 2 ), (x(e) 3 , y(e) 3 ) y (x(e) 4 , y(e) 4 ) están tomados en sentido contrario al movimiento de las manecillas del reloj como se muetra en la figura y el vector b= (b1, b2) es la posición del vector centroide del rectángulo Ωe, ¶ (176) las cuales deberán de ser sustituidas en cada K(e) ij y F (e) j para calcular las integrales en el elemento Ωe. Estas integrales se harán en el programa usando cuadratura Gaussiana, permitiendo reducir el número de cálculos al mínimo pero manteniendo el balance entre precisión y número bajo de operaciones necesarias para realizar las integraciones. Suponiendo que Ω fue dividido en E elementos, estos elementos generan N nodos en total, de los cuales Nd son nodos desconocidos y Nc son nodos conocidos con valor γj , entonces el algoritmo de ensamble de la matriz K y el vector F se puede esquematizar como: Ki,j = (φi, φj) ∀i = 1, 2, ..., E, j = 1, 2, ..., E Fj = (fΩ, φj) ∀j = 1, 2, ..., E ∀j = 1, 2, ..., Nd : bj = bj − γiKi,j ∀i = 1, 2, ..., E Así, se construye una matriz global en la cual están representados los nodos conocidos y los desconocidos, tomando sólo los nodos desconocidos de la matriz K formaremos una matriz A, haciendo lo mismo al vector F formamos el vector b, entonces la solución al problema será la resolución del sistema de ecuaciones lineales Ax= b, este sistema puede resolverse usando por ejemplo el método de gradiente conjugado. El vector x contendrá la solución buscada en los nodos desconocidos Nd. 51 5. Solución de Grandes Sistemas de Ecuaciones En el capítulo anterior se discutió como proceder para transformar un problema de ecuaciones diferenciales parciales con valores en la frontera en un sistema algebraico de ecuaciones y así poder hallar la solución resolviendo el sistema de ecuaciones lineales que se pueden expresar en la forma matricial siguiente Au = b (177) donde la matriz A es bandada (muchos elementos son nulos) y en problemas reales tiene grandes dimensiones. La elección del método específico para resolver el sistema de ecuaciones depende de las propiedades particulares de la matriz A, en las siguientes secciones examinaremos varios métodos y sus implicaciones en cuanto al costo computacional de la resolución del sistema de ecuaciones. En términos generales, si el problema de ecuaciones diferenciales parciales con valores en la frontera en dos dimensiones se discretiza usando una malla de n×m nodos, el sistema algebraico de ecuaciones asociado es del orden de (n × m)2, pero en general la matriz A resultante para el tipo de problemas de interés en el presente trabajo es bandada y definida positiva, por ello es posible hacer uso de estas características para solucionar el sistema algebraico de ecuaciones de forma óptima. Los métodos de resolución del sistema algebraico de ecuaciones Au = b se clasifican en dos grandes grupos: los métodos directos y los métodos iterativos. En los métodos directos la solución u se obtiene en un número fijo de pasos y sólo están sujetos a los errores de redondeo. En los métodos iterativos, se realizan iteraciones para aproximarse a la solución u aprovechando las características propias de la matriz A, tratando de usar un menor número de pasos que en un método directo. Los métodos iterativos rara vez se usan para resolver sistemas lineales de dimensión pequeña (el concepto de dimensión pequeña es muy relativo), ya que el tiempo necesario para conseguir una exactitud satisfactoria rebasa el que requieren los métodos directos. Sin embargo, en el caso de sistemas grandes con un alto porcentaje de elementos cero, son eficientes tanto en el almacenamiento en la computadora como en el tiempo que se invierte en su solución. Por ésta razón al resolver éstos sistemas algebraicos de ecuaciones es preferible aplicar los métodos iterativos tales como son: Jacobi, Gauss-Seidel, sobrerrelajación sucesiva (SOR), etc. Para más información de éstos y otros métodos, así como pruebas en la velocidad de convergencia y precisión, pueden consultarse en [15], [5], [7], [9], [17], [24], [25] y [19]. Cabe hacer mención de que la mayoría del tiempo de cómputo necesario para resolver el problema de ecuaciones diferenciales parciales (EDP), es consumido en la solución del sistema algebraico de ecuaciones asociado a la discretización, por ello es determinante elegir aquel método numérico que minimice el tiempo invertido en este proceso. 52 5.1. Métodos Directos En estos métodos, la solución u se obtiene en un número fijo de pasos y sólo están sujetos a los errores de redondeo. Entre los métodos más importantes podemos encontrar: Eliminación Gausiana, descomposición LU, eliminación bandada y descomposición de Cholesky. Los métodos antes mencionados, se colocaron en orden descendente en cuanto al consumo de recursos computacionales y ascendente en cuanto al aumento en su eficiencia; describiéndose a continuación: Eliminación Gausiana Tal vez es el método más utilizado para encontrar la solución usando métodos directos. Este algoritmo sin embargo no es eficiente, ya que en general, un sistema de N ecuaciones requiere para su almacenaje en memoria de N 2 entradas para la matriz A, pero cerca de N 3/3 + O(N 2) multiplicaciones y N 3/3 + O(N 2) adiciones para encontrar la solución siendo muy costoso computacionalmente. La eliminación Gausiana se basa en la aplicación de operaciones elementales a renglones o columnas de tal forma que es posible obtener matrices equivalentes. En algunos casos nos interesa conocer A−1, por ello si la eliminación se aplica a la matriz aumentada A | I entonces la matriz A de la matriz aumentada se convertirá en la matriz I y la matriz I de la matriz aumentada será A−1. Así, el sistema Au = b se transformará en u = A−1b obteniendo la solución de u. Descomposición LU Sea U una matriz triangular superior obtenida de A por eliminación bandada. Entonces U = L−1A, donde L es una matriz triangular inferior con unos en la diagonal. Las entradas de L−1 pueden obtenerse de los coeficientes mij definidos en el método anterior y pueden ser almacenados estrictamente en las entradas de la diagonal inferior de A ya que estas ya fueron eliminadas. Esto proporciona una factorización LU de A en la misma matriz A ahorrando espacio de memoria. El problema original Au = b se escribe como LU u = b y se reduce a la solución sucesiva de los sistemas lineales triangulares Ly = b y U u = y. (181) La descomposición LU requiere también N 3/3 operaciones aritméticas para la matriz llena, pero sólo N b2 operaciones aritméticas para la matriz con un ancho de banda de b siendo esto más económico computacionalmente. Nótese que para una matriz no singular A, la eliminación de Gausiana (sin redondear filas y columnas) es equivalente a la factorización LU. Eliminación Bandada Cuando se usa la ordenación natural de los nodos, la matriz A que se genera es bandada, por ello se puede ahorrar considerable espacio de almacenamiento en ella. Este algoritmo consiste en triangular a la matriz A por eliminación hacia adelante operando sólo sobre las entradas dentro de la banda central no cero. Así el renglón j es multiplicado por mij = aij /ajj y el resultado es restado al renglón i para i = j + 1, j + 2, .... El resultado es una matriz triangular superior U que tiene ceros abajo de la diagonal en cada columna. Así, es posible resolver el sistema resultante al sustituir en forma inversa las incógnitas. Descomposición de Cholesky Cuando la matriz es simétrica y definida positiva, se obtiene la descomposición LU de la matriz A, así A = LDU = LDLT donde D = diag(U ) es la diagonal con entradas positivas. La mayor ventaja de esta descomposición es que, en el caso en que es aplicable, el costo de cómputo es sustancialmente reducido, ya que requiere de N 3/6 multiplicaciones y N 3/6 adiciones. 5.2. Métodos Iterativos En estos métodos se realizan iteraciones para aproximarse a la solución u aprovechando las características propias de la matriz A, tratando de usar un 54 menor número de pasos que en un método directo, para más información de estos y otros métodos ver [15] y [24]. Un método iterativo en el cual se resuelve el sistema lineal Au = b (182) comienza con una aproximación inicial u0 a la solución u y genera una sucesión de vectores © uk ª ∞k =1 que converge a u. Los métodos iterativos traen consigo un proceso que convierte el sistema Au = b en otro equivalente de la forma u = T u + c para alguna matriz fija T y un vector c. Luego de seleccionar el vector inicial u0 la sucesión de los vectores de la solución aproximada se genera calculando uk = T uk−1 + c ∀k = 1, 2, 3, ... (183) La convergencia a la solución la garantiza el siguiente teorema cuya solución puede verse en [25]. Teorema 27 Si °° T °° < 1, entonces el sistema lineal u = T u+c tiene una solución única u∗ y las iteraciones uk definidas por la fórmula uk = T uk−1+c ∀k = 1, 2, 3, ... convergen hacia la solución exacta u∗ para cualquier aproximación lineal u0. Notemos que mientras menor sea la norma de la matriz T , más rápida es la convergencia, en el caso cuando °° T °° es menor que uno, pero cercano a uno, la convergencia es muy lenta y el número de iteraciones necesario para disminuir el error depende significativamente del error inicial. En este caso, es deseable proponer al vector inicial u0 de forma tal que se mínimo el error inicial. Sin embargo, la elección de dicho vector no tiene importancia si la °° T °° es pequeña ya que la convergencia es rápida. Como es conocido, la velocidad de convergencia de los métodos iterativos dependen de las propiedades espectrales de la matriz de coeficientes del sistema de ecuaciones, cuando el operador diferencial L de la ecuación del problema a resolver es auto-adjunto se obtiene una matriz simétrica y positivo definida y el número de condicionamiento de la matriz A, es por definición cond(A) = λm´ax λm´ın ≥ 1 (184) donde λm´ax y λm´ın es el máximo y mínimo de los eigenvalores de la matriz A. Si el número de condicionamiento es cercano a 1 los métodos numéricos al solucionar el problema convergerá en pocas iteraciones, en caso contrario se requerirán muchas iteraciones. Frecuentemente al usar el método de elemento finito se tiene una velocidad de convergencia de O ¡ 1 h2 ¢ y en el caso de métodos de descomposición de dominio se tiene una velocidad de convergencia de O ¡ 1 h ¢ en el mejor de los casos, donde h es la máxima distancia de separación entre nodos continuos de la partición, es decir, que poseen una pobre velocidad de convergencia cuando h → 0, para más detalles ver [2]. 55 Entre los métodos más usados para el tipo de problemas tratados en el presente trabajo podemos encontrar: Jacobi, Gauss-Seidel, Richardson, relajación sucesiva, gradiente conjugado, gradiente conjugado precondicionado. Los métodos antes mencionados se colocaron en orden descendente en cuanto al consumo de recursos computacionales y ascendente en cuanto al aumento en la eficiencia en su desempeño, describiéndose a continuación: Jacobi Si todos los elementos de la diagonal principal de la matriz A son diferentes de cero aii 6= 0 para i = 1, 2, ...n. Podemos dividir la i−ésima ecuación del sistema lineal (182) por aii para i = 1, 2, ...n, y después trasladamos todas las incógnitas, excepto xi, a la derecha, se obtiene el sistema equivalente u = Bu + d (185) . Las iteraciones del método de Jacobi están definidas por la fórmula donde x(0) i son arbitrarias (i = 1, 2, ....n; k = 1, 2, ....). También el método de Jacobi se puede expresar en términos de matrices. Supongamos por un momento que la matriz A tiene la diagonal unitaria, esto es diag(A) = I. Si descomponemos A = I−B, entonces el sistema dado por la Ecs. (182) se puede reescribir como ¡ I − B ¢ u = b. (187) Para la primera iteración asumimos que k=b; entonces la última ecuación se escribe como u = Bu+k. Tomando una aproximación inicial u0, podemos obtener una mejor aproximación remplazando u por la más resiente aproximación de um. Esta es la idea que subyace en el método Jacobi. El proceso iterativo queda como um+1 = Bum + k. (188) La aplicación del método a la ecuación de la forma Au = b, con la matriz A no cero en los elementos diagonales, se obtiene multiplicando la Ec. (182) por D−1 = £ diag(A) ¤ −1 obteniendo B = I − D−1A, k = D−1b. (189) 56 Gauss-Seidel Este método es una modificación del método Jacobi, en el cual una vez obtenido algún valor de um+1, este es usado para obtener el resto de los valores utilizando los valores más actualizados de um+1. Así, la Ec. (188) Notemos que el método Gauss-Seidel requiere el mismo número de operaciones aritméticas por iteración que el método de Jacobi. Este método se escribe en forma matricial como um+1 = Eum+1 + F um + k (191) donde E y F son las matrices triangular superior e inferior respectivamente. Este método mejora la convergencia con respecto al método de Jacobi en un factor aproximado de 2. Richardson Escribiendo el método de Jacobi como um+1 − um = b − Aum (192) entonces el método Richardson se genera al incorporar la estrategia de sobrerrelajación de la forma siguiente um+1 = um + ω ¡ b − Aum¢ . (193) El método de Richardson se define como um+1 = ¡ I − ωA ¢ um + ωb (194) en la práctica encontrar el valor de ω puede resultar muy costoso computacionalmente y las diversas estrategias para encontrar ω dependen de las características propias del problema, pero este método con un valor ω óptimo resulta mejor que el método de Gauss-Seidel. Relajación Sucesiva Partiendo del método de Gauss-Seidel y sobrerrelajando este esquema, obtenemos y cuando la matriz A es simétrica con entradas en la diagonal positivas, éste método converge si y sólo si A es definida positiva y ω ∈ (0, 2) . En la práctica encontrar el valor de ω puede resultar muy costoso computacionalmente y las diversas estrategias para encontrar ω dependen de las características propias del problema. 57 Gradiente Conjugado El método del gradiente conjugado ha recibido mucha atención en su uso al resolver ecuaciones diferenciales parciales y ha sido ampliamente utilizado en años recientes por la notoria eficiencia al reducir considerablemente en número de iteraciones necesarias para resolver el sistema algebraico de ecuaciones. Aunque los pioneros de este método fueron Hestenes y Stiefel (1952), el interés actual arranca a partir de que Reid (1971) lo planteara como un método iterativo, que es la forma en que se le usa con mayor frecuencia en la actualidad, esta versión está basada en el desarrollo hecho en [9]. La idea básica en que descansa el método del gradiente conjugado consiste en construir una base de vectores ortogonales y utilizarla para realizar la búsqueda de la solución en forma más eficiente. Tal forma de proceder generalmente no sería aconsejable porqué la construcción de una base ortogonal utilizando el procedimiento de Gramm-Schmidt requiere, al seleccionar cada nuevo elemento de la base, asegurar su ortogonalidad con respecto a cada uno de los vectores construidos previamente. La gran ventaja del método de gradiente conjugado radica en que cuando se utiliza este procedimiento, basta con asegurar la ortogonalidad de un nuevo miembro con respecto al último que se ha construido, para que automáticamente esta condición se cumpla con respecto a todos los anteriores. En el algoritmo de gradiente conjugado (CGM), se toman como datos de entrada al sistema Au = b (196) el vector de búsqueda inicial u0 y se calcula r0 = b − Au0, p0 = r0, quedando el método esquemáticamente como: βk+1 = Apk • rk Apk • pk (197) pk+1 = rk − βk+1pk αk+1 = rk • rk Apk+1 • pk+1 uk+1 = uk + αk+1pk+1 rk+1 = rk − αk+1Apk+1. Si denotamos {λi, Vi}N i=1 como las eigensoluciones de A, i.e. AVi = λiVi, i = 1, 2, ..., N. Ya que la matriz A es simétrica, los eigenvalores son reales y podemos ordenarlos por λ1 ≤ λ2 ≤ ... ≤ λN . Definimos el número de condición por Cond(A) = λN /λ1 y la norma de la energía asociada a A por kuk2 A = u •Au 58 El siguiente teorema nos da idea del espectro de convergencia del sistema Au = b para el método de gradiente conjugado. Teorema 28 Sea κ = cond(A) = λm´ax λm´ın ≥ 1, entonces el método de gradiente conjugado satisface la A−norma del error dado por (199) donde em = u − um del sistema Au = b. Notemos que para κ grande se tiene que √κ − 1 √κ + 1 ' 1 − 2 √κ (200) tal que kenkA ' °° e0 °° A exp μ −2 n √κ ¶ (201) de lo anterior podemos esperar un espectro de convergencia del orden de O(√κ) iteraciones, para mayor referencia ver [25]. 5.3. Precondicionadores Una vía que permite mejorar la eficiencia de los métodos iterativos consiste en transformar al sistema de ecuaciones en otro equivalente, en el sentido de que posea la misma solución del sistema original pero que a su vez tenga mejores condiciones espectrales. Esta transformación se conoce como precondicionamiento y consiste en aplicar al sistema de ecuaciones una matriz conocida como precondicionador encargada de realizar el mejoramiento del número de condicionamiento. Una amplia clase de precondicionadores han sido propuestos basados en las características algebraicas de la matriz del sistema de ecuaciones, mientras que por otro lado también existen precondicionadores desarrollados a partir de las características propias del problema que lo origina, un estudio más completo puede encontrarse en [2] y [17]. ¿Qué es un Precondicionador? De una manera formal podemos decir que un precondicionador consiste en construir una matriz C, la cuál es una aproximación en algún sentido de la matriz A del sistema Au = b, de manera tal que si multiplicamos ambos miembros del sistema de ecuaciones original por C−1 obtenemos el siguiente sistema C−1Au = C−1b (202) 59 donde el número de condicionamiento de la matriz del sistema transformado C−1A debe ser menor que el del sistema original, es decir Cond(C−1A) < Cond(A), (203) dicho de otra forma un precondicionador es una inversa aproximada de la matriz original C−1 ' A−1 (204) que en el caso ideal C−1 = A−1 el sistema convergería en una sola iteración, pero el coste computacional del cálculo de A−1 equivaldría a resolver el sistema por un método directo. Se sugiere que C sea una matriz lo más próxima a A sin que su determinación suponga un coste computacional elevado. Dependiendo de la forma de platear el producto de C−1 por la matriz del sistema obtendremos distintas formas de precondicionamiento, estas son: C−1Ax = C−1b Precondicionamiento por la izquierda AC−1Cx = b Precondicionamiento por la derecha C−1 1 AC−1 2 C 2 x = C−1 1 b Precondicionamiento por ambos lados si C puede factorizarse como C = C 1 C 2 . El uso de un precondicionador en un método iterativo provoca que se incurra en un costo de cómputo extra debido a que inicialmente se construye y luego se debe aplicar en cada iteración. Teniéndose que encontrar un balance entre el costo de construcción y aplicación del precondicionador versus la ganancia en velocidad en convergencia del método. Ciertos precondicionadores necesitan poca o ninguna fase de construcción, mientras que otros pueden requerir de un trabajo substancial en esta etapa. Por otra parte la mayoría de los precondicionadores requieren en su aplicación un monto de trabajo proporcional al número de variables; esto implica que se multiplica el trabajo por iteración en un factor constante. De manera resumida un buen precondicionador debe reunir las siguientes características: i) Al aplicar un precondicionador C al sistema original de ecuaciones Au = b, se debe reducir el número de iteraciones necesarias para que la solución aproximada tenga la convergencia a la solución exacta con una exactitud ε prefijada. ii) La matriz C debe ser fácil de calcular, es decir, el costo computacional de la construcción del precondicionador debe ser pequeño comparado con el costo total de resolver el sistema de ecuaciones Au = b. iii) El sistema Cz =r debe ser fácil de resolver. Esto debe interpretarse de dos maneras: a) El monto de operaciones por iteración debido a la aplicación del precondicionador C debe ser pequeño o del mismo orden que las 60 que se requerirían sin precondicionamiento. Esto es importante si se trabaja en máquinas secuenciales. b) El tiempo requerido por iteración debido a la aplicación del precondicionador debe ser pequeño. En computadoras paralelas es importante que la aplicación del precondicionador sea paralelizable, lo cual eleva su eficiencia, pero debe de existir un balance entre la eficacia de un precondicionador en el sentido clásico y su eficiencia en paralelo ya que la mayoría de los precondicionadores tradicionales tienen un componente secuencial grande. El método de gradiente conjugado por si mismo no permite el uso de precondicionadores, pero con una pequeña modificación en el producto interior usado en el método, da origen al método de gradiente conjugado precondicionado que a continuación detallaremos. 5.3.1. Gradiente Conjugado Precondicionado Cuando la matriz A es simétrica y definida positiva se puede escribir como λ1 ≤ uA • u u • u ≤ λn (205) y tomando la matriz C−1 como un precondicionador de A con la condición de que λ1 ≤ uC−1A • u u • u ≤ λn (206) entonces la Ec. (196) se pude escribir como C−1Au = C−1b (207) donde C−1A es también simétrica y definida positiva en el producto interior hu, vi = u • Cv, porque u, C−1Av ® = u • C ¡ C−1Av ¢ (208) = u • Av que por hipótesis es simétrica y definida positiva en ese producto interior. La elección del producto interior h•, •i quedará definido como hu, vi = u • C−1Av (209) por ello las Ecs. (197[1]) y (197[3]), se convierten en αk+1 = rk • rk pk+1 • C−1pk+1 (210) 61 y βk+1 = pk • C−1rk pk • Apk (211) generando el método de gradiente conjugado precondicionado con precondicionador C−1. Es necesario hacer notar que los métodos gradiente conjugado y gradiente conjugado precondicionado sólo difieren en la elección del producto interior. Para el método de gradiente conjugado precondicionado, los datos de entrada son un vector de búsqueda inicial u0 y el precondicionador C−1. Calculándose r0 = b − Au0, p = C−1r0, quedando el método esquemáticamente como: βk+1 = pk • C−1rk pk • Apk (212) pk+1 = rk − βk+1pk αk+1 = rk • rk pk+1 • C−1pk+1 uk+1 = uk + αk+1pk+1 rk+1 = C−1rk − αk+1Apk+1. Algoritmo Computacional del Método Dado el sistema Au = b, con la matriz A simétrica y definida positiva de dimensión n×n. La entrada al método será una elección de u0 como condición inicial, ε > 0 como la tolerancia del método, N como el número máximo de iteraciones y la matriz de precondicionamiento C−1 de dimensión n × n, el algoritmo del método de gradiente La salida del método será la solución aproximada x = (x1, ..., xN ) y el residual r = (r1, ..., rN ). En el caso del método sin precondicionamiento, C−1 es la matriz identidad, que para propósitos de optimización sólo es necesario hacer la asignación de vectores correspondiente en lugar del producto de la matriz por el vector. En el caso de que la matriz A no sea simétrica, el método de gradiente conjugado puede extenderse para soportarlas, para más información sobre pruebas de convergencia, resultados numéricos entre los distintos métodos de solución del sistema algebraico Au = b generada por la discretización de un problema elíptico y como extender estos para matrices no simétricas ver [9] y [7]. Clasificación de los Precondicionadores En general se pueden clasificar en dos grandes grupos según su manera de construcción: los algebraicos o a posteriori y los a priori o directamente relacionados con el problema continuo que lo origina. 5.3.2. Precondicionador a Posteriori Los precondicionadores algebraicos o a posteriori son los más generales, ya que sólo dependen de la estructura algebraica de la matriz A, esto quiere decir que no tienen en cuenta los detalles del proceso usado para construir el sistema de ecuaciones lineales Au = b. Entre estos podemos citar los métodos de precondicionamiento del tipo Jacobi, SSOR, factorización incompleta, inversa aproximada, diagonal óptimo y polinomial. Precondicionador Jacobi El método precondicionador Jacobi es el precondicionador más simple que existe y consiste en tomar en calidad de precondicionador a los elementos de la diagonal de A Cij = ⎧⎨ ⎩ Aij si i = j 0 si i 6= j. (213) Debido a que las operaciones de división son usualmente más costosas en tiempo de cómputo, en la práctica se almacenan los recíprocos de la diagonal de A. Ventajas: No necesita trabajo para su construcción y puede mejorar la convergencia. Desventajas: En problemas con número de condicionamiento muy grande, no es notoria la mejoría en el número de iteraciones. 63 Precondicionador SSOR Si la matriz original es simétrica, se puede descomponer como en el método de sobrerrelajamiento sucesivo simétrico (SSOR) de la siguiente manera A = D + L + LT (214) donde D es la matriz de la diagonal principal y L es la matriz triangular inferior. (215) en la práctica la información espectral necesaria para hallar el valor óptimo de ω es demasiado costoso para ser calculado. Ventajas: No necesita trabajo para su construcción, puede mejorar la convergencia significativamente. Desventajas: Su paralelización depende fuertemente del ordenamiento de las variables. Precondicionador de Factorización Incompleta Existen una amplia clase de precondicionadores basados en factorizaciones incompletas. La idea consiste en que durante el proceso de factorización se ignoran ciertos elementos diferentes de cero correspondientes a posiciones de la matriz original que son nulos. La matriz precondicionadora se expresa como C = LU , donde L es la matriz triangular inferior y U la superior. La eficacia del método depende de cuán buena sea la aproximación de C−1 con respecto a A−1. El tipo más común de factorización incompleta se basa en seleccionar un subconjunto S de las posiciones de los elementos de la matriz y durante el proceso de factorización considerar a cualquier posición fuera de éste igual a cero. Usualmente se toma como S al conjunto de todas las posiciones (i, j) para las que Aij 6= 0. Este tipo de factorización es conocido como factorización incompleta LU de nivel cero, ILU(0). El proceso de factorización incompleta puede ser descrito formalmente como sigue: Para cada k, si i, j > k: Sij = ⎧⎨ ⎩ Aij − AijA−1 ij Akj Si (i, j) ∈ S Aij Si (i, j) / ∈ S. (216) Una variante de la idea básica de las factorizaciones incompletas lo constituye la factorización incompleta modificada que consiste en que si el producto Aij − AijA−1 ij Akj 6= 0 (217) y el llenado no está permitido en la posición (i, j), en lugar de simplemente descartarlo, esta cantidad se le substrae al elemento de la diagonal Aij .Matemáticamente esto corresponde a forzar a la matriz precondicionadora a tener la misma suma por filas que la matriz original. Esta variante resulta de interés puesto 64 que se ha probado que para ciertos casos la aplicación de la factorización incompleta modificada combinada con pequeñas perturbaciones hace que el número de condicionamiento espectral del sistema precondicionado sea de un orden inferior. Ventaja: Puede mejorar el condicionamiento y la convergencia significativamente. Desventaja: El proceso de factorización es costoso y difícil de paralelizar en general. Precondicionador de Inversa Aproximada El uso del precondicionador de inversas aproximada se ha convertido en una buena alternativa para los precondicionadores implícitos debido a su naturaleza paralelizable. Aquí se construye una matriz inversa aproximada usando el producto escalar de Frobenius. Sea S ⊂ Cn, el subespacio de las matrices C donde se busca una inversa aproximada explícita con un patrón de dispersión desconocido. La formulación del problema esta dada como: Encontrar C 0 ∈ S tal que C 0 = argm´ınC∈S °° AC − I °° . (218) Además, esta matriz inicial C 0 puede ser una inversa aproximada de A en un sentido estricto, es decir, °° ° AC 0 − I °°° = ε < 1. (219) Existen dos razones para esto, primero, la ecuación (219) permite asegurar que C 0 no es singular (lema de Banach), y segundo, esta será la base para construir un algoritmo explícito para mejorar C 0 y resolver la ecuación Au = b. La construcción de C 0 se realiza en paralelo, independizando el cálculo de cada columna. El algoritmo permite comenzar desde cualquier entrada de la columna k, se acepta comúnmente el uso de la diagonal como primera aproximación. Sea rk el residuo correspondiente a la columna k-ésima, es decir rk = ACk − ek (220) y sea Ik el conjunto de índices de las entradas no nulas en rk , es decir, Ik = {i = {1, 2, ..., n} | rik 6= 0} . Si Lk = {l = {1, 2, ..., n} | Clk 6= 0} , entonces la nueva entrada se busca en el conjunto Jk = {j ∈ Lc k | Aij 6= 0, ∀i ∈ Ik} . En realidad las únicas entradas consideradas en Ck son aquellas que afectan las entradas no nulas de rk Esta estrategia define el nuevo índice seleccionado jk atendiendo solamente al conjunto Lk , lo que nos lleva a un nuevo óptimo donde se actualizan todas las entradas correspondientes a los índices de Lk . Esto mejora el criterio de (218) donde el nuevo índice se selecciona manteniendo las entradas correspondientes a los índices de Lk . donde ˜ Cl es el vector con entradas no nulas ik h (1 ≤ h ≤ l) . Cada una de ellas se obtiene evaluado el determinante correspondiente que resulta de remplazar la última fila del det ³ Gk l ´ por et h, con 1 ≤ l ≤ pk . Evidentemente, los cálculos de °° ACk − ek ° ° 2 2 y de Ck pueden actualizarse añadiendo la contribución de la última entrada j ∈ Jk a la suma previa de 1 a pk −1. En la práctica, det ³ Gk l ´ se calcula usando la descomposición de Cholesky puesto que Gk l es una matriz simétrica y definida positiva. Esto sólo involucra la factorización de la última fila y columna si aprovechamos la descomposición de Gk l−1 necesitándose solamente una sustitución por descenso. Finalmente, para obtener ˜ Cl debe resolverse el sistema Gk l vl = el, con ˜ Cik1 l = vhl, (1 ≤ h ≤ l) . Ventaja: Puede mejorar el condicionamiento y la convergencia significativamente y es fácilmente paralelizable. Desventaja: El proceso construcción es algo laborioso. 5.3.3. Precondicionador a Priori Los precondicionadores a priori son más particulares y dependen para su construcción del conocimiento del proceso de discretización de la ecuación diferencial parcial, dicho de otro modo dependen más del proceso de construcción de la matriz A que de la estructura de la misma. Estos precondicionadores usualmente requieren de más trabajo que los del tipo algebraico discutidos anteriormente, sin embargo permiten el desarrollo de métodos de solución especializados más rápidos que los primeros. 66 Veremos algunos de los métodos más usados relacionados con la solución de ecuaciones diferenciales parciales en general y luego nos concentraremos en el caso de los métodos relacionados directamente con descomposición de dominio. En estos casos el precondicionador C no necesariamente toma la forma simple de una matriz, sino que debe ser visto como un operador en general. De aquí que C podría representar al operador correspondiente a una versión simplificada del problema con valores en la frontera que deseamos resolver. Por ejemplo se podría emplear en calidad de precondicionador al operador original del problema con coeficientes variables tomado con coeficientes constantes. En el caso del operador de Laplace se podría tomar como precondicionador a su discretización en diferencias finitas centrales. Por lo general estos métodos alcanzan una mayor eficiencia y una convergencia óptima, es decir, para ese problema en particular el precondicionador encontrado será el mejor precondicionador existente, llegando a disminuir el número de iteraciones hasta en un orden de magnitud. Donde muchos de ellos pueden ser paralelizados de forma efectiva. El Uso de la Parte Simétrica como Precondicionador La aplicación del método del gradiente conjugado en sistemas no auto-adjuntos requiere del almacenamiento de los vectores previamente calculados. Si se usa como precondicionador la parte simétrica (A + AT )/2 (223) de la matriz de coeficientes A, entonces no se requiere de éste almacenamiento extra en algunos casos, resolver el sistema de la parte simétrica de la matriz A puede resultar más complicado que resolver el sistema completo. El Uso de Métodos Directos Rápidos como Precondicionadores En muchas aplicaciones la matriz de coeficientes A es simétrica y positivo definida, debido a que proviene de un operador diferencial auto-adjunto y acotado. Esto implica que se cumple la siguiente relación para cualquier matriz B obtenida de una ecuación diferencial similar c1 ≤ xT Ax xT Bx ≤ c2 ∀x (224) donde c1 y c2 no dependen del tamaño de la matriz. La importancia de esta propiedad es que del uso de B como precondicionador resulta un método iterativo cuyo número de iteraciones no depende del tamaño de la matriz. La elección más común para construir el precondicionador B es a partir de la ecuación diferencial parcial separable. El sistema resultante con la matriz B puede ser resuelto usando uno de los métodos directos de solución rápida, como pueden ser por ejemplo los basados en la transformada rápida de Fourier. Como una ilustración simple del presente caso obtenemos que cualquier operador elíptico puede ser precondicionado con el operador de Poisson. 67 Construcción de Precondicionadores para Problemas Elípticos Empleando DDM Existen una amplia gama de este tipo de precondicionadores, pero son específicos al método de descomposición de dominio usado, para el método de subestructuración, los más importantes se derivan de la matriz de rigidez y por el método de proyecciones, el primero se detalla en la sección (6.2.1) y el segundo, conjuntamente con otros precondicionadores pueden ser consultados en [11], [5], [4] y [2]. La gran ventaja de este tipo de precondicionadores es que pueden ser óptimos, es decir, para ese problema en particular el precondicionador encontrado será el mejor precondicionador existente, llegando a disminuir el número de iteraciones hasta en un orden de magnitud. 68 6. Métodos de Descomposición de Dominio (DDM) La solución numérica por los esquemas tradicionales de discretización tipo elemento finito y diferencias finitas generan una discretización del problema, la cual es usada para generar un sistema de ecuaciones algebraicos Au = b. Este sistema algebraico en general es de gran tamaño para problemas reales, al ser estos algoritmos secuenciales su implantación suele hacerse en equipos secuenciales y por ello no es posible resolver muchos problemas que involucren el uso de una gran cantidad de memoria, actualmente para tratar de subsanar dicha limitante, se usa equipo paralelo para soportar algoritmos secuenciales, haciendo ineficiente su implantación en dichos equipos. Los métodos de descomposición de dominio son un paradigma natural usado por la comunidad de modeladores. Los sistemas físicos son descompuestos en dos o más subdominios contiguos basados en consideraciones fenomenológicas. Estas descomposiciones basadas en dominios físicos son reflejadas en la ingeniería de software del código correspondiente. Los métodos de descomposición de dominio permiten tratar los problemas de tamaño considerable, empleando algoritmos paralelos en computadoras secuenciales y/o paralelas. Esto es posible ya que cualquier método de descomposición de dominio se basa en la suposición de que dado un dominio computacional Ω, este se puede particionar en subdominios Ωi,i = 1, 2, ..., E entre los cuales puede o no existir traslape. Entonces el problema es reformulado en términos de cada subdominio (empleando algún método del tipo elemento finito) obteniendo una familia de subproblemas de tamaño reducido independientes en principio entre si, que están acoplados a través de la solución en la interfaz de los subdominios que es desconocida. De esta manera, podemos clasificar de manera burda a los métodos de descomposición de dominio, como aquellos en que: existe traslape entre los subdominios y en los que no existe traslape. A la primera clase pertenece el método de Schwarz (en el cual el tamaño del traslape es importante en la convergencia del método) y a los de la segunda clase pertenecen los métodos del tipo subestructuración (en el cual los subdominios sólo tienen en común los nodos de la frontera interior). La computación en paralelo es una técnica que nos permite distribuir una gran carga computacional entre muchos procesadores. Y es bien sabido que una de las mayores dificultades del procesamiento en paralelo es la coordinación de las actividades de los diferentes procesadores y el intercambio de información entre los mismos [21] mediante el paso de mensajes. Así, mediante los métodos de descomposición de dominio, la programación orientada a objetos y esquemas de paralelización que usan el paso de mensajes, es posible construir aplicaciones que coadyuven a la solución de problemas concomitantes en ciencia e ingeniería, ya que permiten utilizar todas las capacidades del cómputo en paralelo (supercomputadoras, clusters o grids), de esta forma es posible atacar una gran variedad de problemas que sin estas técnicas es imposible hacerlo de manera flexible y eficiente. Pero hay que notar que existen una amplia gama de problemas que nos 69 interesan resolver que superan las capacidades de cómputo actuales, ya sea por el tiempo requerido para su solución, por el consumo excesivo de memoria o ambos. La lista de los métodos de descomposición de dominio y el tipo de problemas que pueden ser atacados por estos, es grande y está en constante evolución, ya que se trata de encontrar un equilibrio entre la complejidad del método (aunada a la propia complejidad del modelo), la eficiencia en el consumo de los recursos computacionales y la precisión esperada en la solución encontrada por los diversos métodos y las arquitecturas paralelas en la que se implante. A continuación describiremos algunos de estos métodos generales. En este capítulo se considerarán problemas con valor en la frontera (VBVP) de la forma Lu = f en Ω (225) u = g en ∂Ω donde Lu = −∇ • a • ∇u + cu (226) como un caso particular del operador elíptico definido por la Ec. (79) de orden dos. 6.1. Método de Schwarz El método fue desarrollado por Hermann Amandus Schwarz en 1869 (no como un método de descomposición de dominio), ya que en esos tiempos los matemáticos podían resolver problemas con geometrías sencillas de manera analítica, pero no tenían una idea clara de como poder resolver problemas que involucraran el traslape de esas geometrías sencillas. Como se conocía la solución para las geometrías sencillas por separado, la idea de Schwarz fue usar estas para conocer la solución en la geometría resultante al tener traslape, para más detalle ver [1]. Para describir el método, consideremos primero un dominio Ω que está formado de dos subdominios Ω1 y Ω2 traslapados, es decir Ω1 ∩ Ω2 6= ∅, entonces Ω = Ω1 ∪ Ω2 y denotemos a Σ1 = ∂Ω1 ∩ Ω2, Σ2 = ∂Ω2 ∩ Ω1 y Ω1,2 = Ω1 ∩ Ω2, como se muestra en la figura para dos dominios distintos: Figura 3: Dominio Ω subdividido en dos subdominios Ω1 y Ω2. 70 La forma original del método iterativo de Schwarz conocido como métodos alternantes de Schwarz, consiste en resolver sucesivamente los siguientes problemas. Sea uo una función de inicialización definida en Ω, que se nulifica en ∂Ω, (228) resolviendo los problemas secuencialmente en cada subdominio (por ejemplo con el método de elemento finito). Este método se conoce como Schwarz multiplicativo. El método alternante de Schwarz dado por las Ecs. (227) y (228) converge a la solución u de (225) si suponemos alguna suavidad en los subdominios Ω1 y Ω2, ya que existen constantes C1 y C2 ∈ (0, 1) tal que para todo k ≥ 0 se tiene °° las constantes C1 y C2 de reducción de error deben de estar bastante cerca de 1 si la región de traslape Ω1,2 es delgada, la prueba de esta estimación puede encontrarse en [12]. (231) resolviendo los problemas en paralelo de cada subdominio (por ejemplo con el método de elemento finito). Este método se conoce como Schwarz aditivo. La convergencia de este método en general requiere de algunas hipótesis adicionales, pero si converge, el número de iteraciones necesarias para converger será del doble que el método Schwarz multiplicativo. 71 La generalización del método de Schwarz en el caso en que Ω es particionada en E > 2 subdominios traslapados puede describirse como sigue: Descomponiendo el dominio Ω en E subdominios Ωe con traslape como por ejemplo, la descomposición siguiente Figura 4: Descomposición de Ω en múltiples subdominios con traslape para el método de Schwarz. Entonces para resolver el problema por el método de Schwarz, primeramente es necesario definir los siguientes subespacios del espacio de Sobolev H1(Ω), Observaciones: • La precisión del método depende fuertemente del número de iteraciones realizadas en el proceso iterativo y converge a la precisión usada en la solución de cada subdominio en el mejor de los casos. • El método aditivo de Schwarz es secuencial, en el caso del método multiplicativo de Schwarz es paralelizable pero tiene una parte serial importante en el algoritmo y su convergencia no es la óptima en esta formulación, pero existen variantes del método que permiten remediar esto, para más detalles ver [2], [4] y [5]. • Hay que notar que por cada subdominio (supóngase n) y en cada iteración (supóngase I) se resuelve un problema para cada Ωi, esto significa que si se usa el método de elemento finito para resolver el problema local donde se usen en promedio r iteraciones para resolver el sistema lineal (no tomando en cuenta el costo invertido en generar las matrices), el total de iteraciones necesarias para resolver el problema en el dominio Ω será r ∗ n ∗ I, resultando muy costoso computacionalmente con respecto a otros métodos de descomposición de dominio. 73 6.2. Método de Subestructuración Consideremos el problema dado por la Ec. (225) en el dominio Ω, el cual es subdividido en E subdominios Ωi, i = 1, 2, ..., E sin traslape, es decir Ωi ∩ Ωj = ∅ ∀i 6= j y Ω = [E i=1 Ωi, (241) y al conjunto Σ = [E i=1 Σi, si Σi = ∂Ωi\∂Ω (242) lo llamaremos la frontera interior de los subdominios, un ejemplo se muestra en la figura: Figura 5: Dominio Ω descompuesto en subdominios Ωi, con i = 1, 2, ..., 9. Sin perdida de generalidad tomemos g = 0 en ∂Ω, notemos que siempre es posible poner el problema de la Ec. (225) como uno con condiciones de frontera Dirichlet que se nulifiquen mediante la adecuada manipulación del término del lado derecho de la ecuación. Primeramente sea D ⊂ H1 0 (Ω) un espacio lineal de funciones de dimensión finita N, en el cual esté definido un producto interior denotado para cada u, v ∈ D por u • v = hu, vi (243) además sean ˜DI , ˜DΣ y ¯DΣ subespacios lineales de D con la propiedad de que ˜D Σ y ¯DΣ cada uno por separado sean los complementos algebraicos con respecto a D de ˜DI . Más explícitamente D = ˜DI + ˜DΣ y ˜DI ∩ ˜DΣ = {0} (244) y D = ˜DI + ¯DΣ y ˜DI ∩ ¯DΣ = {0} . (245) 74 Adicionalmente ¯DΣ es ortogonal a ˜DI , i.e. ¯DΣ es el complemento ortogonal de ˜DI , como se muestra en la figura: Figura 6: Esquematización de los subespacios ˜DI , ˜DΣ y ¯DΣ del espacio D. Sean ξI = © wiI | i = 1, 2, ..., NI ª , ξΣ = {wαΣ | α = 1, 2, ..., NΣ} (246) bases linealmente independientes de ˜DI y ˜DΣ respectivamente, tales que ξ = ξI ∪ ξΣ (247) sea una base del espacio D, donde N = NI + NΣ y sin perdida de generalidad podemos suponer que ξ es una base ordenada, es decir, primero están los vectores linealmente independientes de ξI y después los de ξΣ. Denotaremos al dual algebraico de D por D∗, el cual será el espacio lineal de funciones definidas en D. Definimos también los operadores proyección P1, P2 y P3 de ˜DI , ¯DΣ y ˜DΣ respectivamente, notemos que podemos definir una biyección entre los subespacios ¯D Σ y ˜DΣ por medio de las proyecciones P2 : ˜DΣ → ¯DΣ y P3 : ¯DΣ → ˜DΣ, al tener la misma cardinalidad y ser complementos algebraicos con respecto a D de ˜DI . Como trabajaremos con espacios de dimensión finita, cuando ξI ⊂ ˜DI y ξΣ ⊂ ˜DΣ se mantengan fijos, se puede asociar con cada u ∈ D un único vector u≡ ______(uI , uΣ) donde uI = ³ uI1 , ..., uINI ´ y uΣ = ³ uΣ1 , ..., uΣNΣ ´ , entonces u= (u1, ..., uN ) ∈ RN tal que u = XN i=1 uiwi (248) donde wi ∈ ξ y ui ∈ R para i = 1, ..., N. En especial cuando uΣ ∈ ˜DΣ y uI ∈ ˜DI 75 podemos asociar un único vector uΣ∈ RNΣ y uI ∈ RNI respectivamente tal que uΣ = XN α=1 uαwαΣ y uI = XN i=1 uiwiI . (249) Las funciones de D → RN definidas de está manera son una biyeción a la cual nos referiremos como la inmersión natural de RN en D. Adicionalmente, las imágenes bajo la inmersión natural de ˜DI y ˜DΣ son isomorfos a RNI y RNΣ respectivamente, sin perdida de generalidad podemos suponer que la base (u1, ..., uN ) ∈ RN es una base ordenada, es decir, suponemos que primero aparecen los vectores de (v1, ..., vNI ) ∈ RNI y después los vectores (w1, ..., wNΣ ) ∈ RNΣ , entonces la base para RN será (v1, ..., vNI , w1, ..., wNΣ ) . El producto interior Euclidiano en RN , para cada par u∈ RN y v∈ RN , denotado por u•v, será definido por u • v ≡ XN i=1 uivi, (250) el cual no necesariamente coincide con el producto interior del espacio D definido en la Ec. (243). Si adicionalmente suponemos que existe una familia ortogonal ˜DIi (Ωi) de subespacios linealmente independientes del subespacio de ˜DI (Ω), con i = 1, ..., E, es decir n ˜D I1 (Ω1) , ..., ˜DIE (ΩE) o tales que ˜D I = XE i=1 ˜D Ii (251) y denotamos por P1i , i = 1, ..., E, al operador proyección sobre ˜DIi , entonces P1 = XE i=1 P1i (252) además, para cada j = 1, ..., E, sea ξIj ≡ n wj I | i = 1, ..., E o (253) tal que ξIj es una base de ˜DIj (las funciones wj I pueden ser las φi usadas en el método de elemento finito o cualquier otro tipo de funciones base es construida sumando las AΣΣ i según el orden de los nodos globales versus los nodos locales. podrían ser construida sumando las Si y bi respectivamente según el orden de los nodos globales versus los nodos locales. 78 El sistema lineal virtual obtenido de esta forma (268) se resuelve eficientemente usando el método de gradiente conjugado visto en la sección (5.2), para ello no es necesario construir la matriz S con las contribuciones de cada Si correspondientes al subdominio i, lo que hacemos es pasar a cada subdominio el vector uΣ i correspondiente a la i-ésima iteración del método de gradiente conjugado para que en cada subdominio se evalué ˜uΣ i = S i uΣ i localmente y con el resultado se forma el vector ˜uΣ = PE i=1 ˜uΣ i y se continué con los demás pasos del método. Esto es ideal para una implementación en paralelo del método de gradiente conjugado. Observación 29 Notemos que el normalmente las matrices locales S i y ³ AII i ´ −1 no se construyen, ya que estas serian matrices densas y su construcción es computacionalmente muy costosa. Y como sólo nos interesa el producto SyΣ, o más precisamente hPE i=1 S i i yΣ, entonces si llamamos yΣi al vector correspondiente al subdominio i, entonces tendremos z = ³ AΣΣ − AΣI ¡ AII¢−1 AIΣ ´ yΣi . (270) Para evaluar eficientemente esta expresión, realizamos las siguientes operaciones equivalentes x1 = AΣΣyΣi (271) x2 = ³ AΣI ¡ AII¢−1 AIΣ ´ yΣi z = x1 − x2 la primera y tercera expresión no tienen ningún problema en su evaluación, para la segunda expresión tendremos que hacer x3 = AIΣyΣi (272) con este resultado intermedio, deberíamos calcular x4 = ¡ AII¢−1 x3 (273) pero como no contamos con ¡ AII¢−1 , entonces multiplicamos la expresión por AII obteniendo AII x4 = AII ¡ AII¢−1 x3 (274) al simplificar, tenemos AII x4 = x3. (275) Esta última expresión puede ser resuelta usando Factorización LU o Gradiente Conjugado (cada una de estas opciones tiene ventajas y desventajas que deben ser evaluadas al momento de implementar el código para un problema particular). Una vez obtenido x3, podremos calcular x2 = AΣI x3 (276) completando la secuencia de operaciones necesaria para obtener S i yΣ. 79 Una vez resuelto el sistema de la Ec. (269) en el que hemos encontrado la solución para los nodos de la frontera interior uΣ, entonces debemos resolver localmente los uI i correspondientes a los nodos interiores para cada subespacio Ωi, para esto empleamos uI i = ³ AII i ´ −1 ³ bI i − AIΣ i uΣi ´ (277) para cada i = 1, 2, ..., E, quedando así resuelto el problema Au = b tanto en los nodos interiores uI i como en los de la frontera interior uΣi correspondientes a cada subespacio Ωi. Observación 30 En la evaluación de uI i = ³ AII i ´ −1 ³ bI i − AIΣ i uΣi ´ , esta nuevamente involucrado ³ AII i ´ −1 , por ello deberemos de usar el siguiente procedimiento para evaluar eficientemente esta expresión, realizando las siguientes operaciones equivalentes x4 = bI i − AIΣ i uΣi (278) uI i = ³ AII i ´ −1 x4 multiplicando por AII i a la última expresión, obtenemos AII i uI i = AII i ³ AII i ´ −1 x4 (279) simplificando, tenemos AII i uI i = x4 (280) esta última expresión puede ser resuelta usando Factorización LU o Gradiente Conjugado. Como se indico en las dos últimas observaciones, para resolver el sistema AII i x = b podemos usar Factorización LU, Gradiente Conjugado o cualquier otro método para resolver sistemas lineales, pero deberá de usarse aquel que proporcione la mayor velocidad en el cálculo o que consuma la menor cantidad de memoria (ambas condicionantes son mutuamente excluyentes), por ello la decisión de que método usar deberá de tomarse al momento de tener que resolver un problema particular en un equipo dado. Para usar el método de Factorización LU, se deberá primeramente de factorizar la matriz bandada AII i en una matriz LU que es densa (del orden de (Nδ)2 por ello consume una gran cantidad de memoria) pero esta operación sólo se deberá de realizar una vez por cada subdominio, y para solucionar los diversos sistemas lineales AII i x = b sólo será necesario evaluar los sistemas Ly = b (281) U x = y 80 en donde y es un vector auxiliar. Esto proporciona una manera muy eficiente de evaluar el sistema lineal pero el consumo en memoria para un problema particular puede ser excesivo. Por ello si el problema involucra una gran cantidad de nodos interiores y el equipo en el que se implantará la ejecución del programa tiene una cantidad de memoria muy limitada, es recomendable usar el método de Gradiente Conjugado, este consume una cantidad de memoria adicional muy pequeña y pese a que no es tan eficiente (dos o tres veces mas lento) como la Factorización LU, si proporciona un buen desempeño versus el obtenido al realizar el cálculo directo de ³ AII i ´ −1 . De esta forma, es posible adaptar el código para tomar en cuenta la implementación de este en un equipo de cómputo en particular y poder sacar el máximo provecho al método de Subestructuración en la resolución de problemas elípticos de gran envergadura. En lo que resta del presente trabajo, se asume que el método empleado para resolver AII i x = b en sus respectivas variantes necesarias para evitar el cálculo de ³ AII i ´ −1 es la Factorización LU, logrando así el máximo desempeño en velocidad en tiempo de ejecución. Para más detalles de la forma como se presento este método ver [10] y para ver otras formulaciones ver [5], [4] y [2]. Observaciones: • La precisión del método es la misma que la usada por los interpoladores en la descomposición de los subdominios. • El método es paralelizable permitiendo que cada subdominio sea manipulado por un procesador y además es posible usar en cada subdominio para el manejo de las matrices generadas diversos esquemas de paralelización para aumentar el rendimiento de los procesadores. • Hay que notar que por cada subdominio (supóngase E) se resuelven sólo los nodos de la frontera interior k, esto significa que en promedio se usan menos de k iteraciones en el método de gradiente conjugado para resolver el sistema lineal asociado (estas iteraciones son en cantidad menor que las realizadas por el método de Schwarz para un sólo subdominio, si Ωi es tomado aproximadamente del mismo tamaño en ambos métodos). Y como la solución de los nodos interiores no conllevan iteraciones adicionales, este método resulta más económico computacionalmente que el método de Schwarz. Está economía se hace aún más patente cuando se usa un buen precondicionador, logrando bajar hasta en un orden de magnitud el número de iteraciones del caso no precondicionado. 81 6.2.1. Precondicionador Derivado de la Matriz de Rigidez En esta sección detallaremos la construcción del precondicionador derivado de la matriz de rigidez para problemas elípticos usando en el método de subestructuración. 7. El Cómputo en Paralelo Los sistemas de cómputo con procesamiento en paralelo surgen de la necesidad de resolver problemas complejos en un tiempo razonable, utilizando las ventajas de memoria, velocidad de los procesadores, formas de interconexión de estos y distribución de la tarea, a los que en su conjunto denominamos arquitectura en paralelo. Entenderemos por una arquitectura en paralelo a un conjunto de procesadores interconectados capaces de cooperar en la solución de un problema. Así, para resolver un problema en particular, se usa una o combinación de múltiples arquitecturas (topologías), ya que cada una ofrece ventajas y desventajas que tienen que ser sopesadas antes de implementar la solución del problema en una arquitectura en particular. También es necesario conocer los problemas a los que se enfrenta un desarrollador de programas que se desean correr en paralelo, como son: el partir eficientemente un problema en múltiples tareas y como distribuir estas según la arquitectura en particular con que se trabaje. 7.1. Arquitecturas de Software y Hardware En esta sección se explican en detalle las dos clasificaciones de computadoras más conocidas en la actualidad. La primera clasificación, es la clasificación clásica de Flynn en dónde se tienen en cuenta sistemas con uno o varios procesadores, la segunda clasificación es moderna en la que sólo tienen en cuenta los sistemas con más de un procesador. El objetivo de esta sección es presentar de una forma clara los tipos de clasificación que existen en la actualidad desde el punto de vista de distintos autores, así como cuáles son las ventajas e inconvenientes que cada uno ostenta, ya que es común que al resolver un problema particular se usen una o más arquitecturas de hardware interconectadas generalmente por red. 7.1.1. Clasificación de Flynn Clasificación clásica de arquitecturas de computadoras que hace alusión a sistemas con uno o varios procesadores, Flynn la publicó por primera vez en 1966 y por segunda vez en 1970. Esta taxonomía se basa en el flujo que siguen los datos dentro de la máquina y de las instrucciones sobre esos datos. Se define como flujo de instrucciones al conjunto de instrucciones secuenciales que son ejecutadas por un único procesador y como flujo de datos al flujo secuencial de datos requeridos por el flujo de instrucciones. Con estas consideraciones, Flynn clasifica los sistemas en cuatro categorías: Single Instruction stream, Single Data stream (SISD) Los sistemas de este tipo se caracterizan por tener un único flujo de instrucciones sobre un único flujo de datos, es decir, se ejecuta una instrucción detrás de otra. Este es el concepto de arquitectura serie de Von Neumann donde, en cualquier 87 momento, sólo se ejecuta una única instrucción, un ejemplo de estos sistemas son las máquinas secuenciales convencionales. Single Instruction stream, Multiple Data stream (SIMD) Estos sistemas tienen un único flujo de instrucciones que operan sobre múltiples flujos de datos. Ejemplos de estos sistemas los tenemos en las máquinas vectoriales con hardware escalar y vectorial. El procesamiento es síncrono, la ejecución de las instrucciones sigue siendo secuencial como en el caso anterior, todos los elementos realizan una misma instrucción pero sobre una gran cantidad de datos. Por este motivo existirá concurrencia de operación, es decir, esta clasificación es el origen de la máquina paralela. El funcionamiento de este tipo de sistemas es el siguiente. La unidad de control manda una misma instrucción a todas las unidades de proceso (ALUs). Las unidades de proceso operan sobre datos diferentes pero con la misma instrucción recibida. Existen dos alternativas distintas que aparecen después de realizarse esta clasificación: Arquitectura Vectorial con segmentación, una CPU única particionada en unidades funcionales independientes trabajando sobre flujos de datos concretos. Arquitectura Matricial (matriz de procesadores), varias ALUs idénticas a las que el procesador da instrucciones, asigna una única instrucción pero trabajando sobre diferentes partes del programa. Multiple Instruction stream, Single Data stream (MISD) Sistemas con múltiples instrucciones que operan sobre un único flujo de datos. Este tipo de sistemas no ha tenido implementación hasta hace poco tiempo. Los sistemas MISD se contemplan de dos maneras distintas: Varias instrucciones operando simultáneamente sobre un único dato. Varias instrucciones operando sobre un dato que se va convirtiendo en un resultado que será la entrada para la siguiente etapa. Se trabaja de forma segmentada, todas las unidades de proceso pueden trabajar de forma concurrente. Multiple Instruction stream,Multiple Data stream (MIMD) Sistemas con un flujo de múltiples instrucciones que operan sobre múltiples datos. Estos sistemas empezaron a utilizarse antes de la década de los 80s. Son sistemas con memoria compartida que permiten ejecutar varios procesos simultáneamente (sistema multiprocesador). Cuando las unidades de proceso reciben datos de una memoria no compartida estos sistemas reciben el nombre de MULTIPLE SISD (MSISD). En 88 arquitecturas con varias unidades de control (MISD Y MIMD), existe otro nivel superior con una unidad de control que se encarga de controlar todas las unidades de control del sistema (ejemplo de estos sistemas son las máquinas paralelas actuales). 7.1.2. Categorías de Computadoras Paralelas Clasificación moderna que hace alusión única y exclusivamente a los sistemas que tienen más de un procesador (i.e máquinas paralelas). Existen dos tipos de sistemas teniendo en cuenta su acoplamiento: Los sistemas fuertemente acoplados son aquellos en los que los procesadores dependen unos de otros. Los sistemas débilmente acoplados son aquellos en los que existe poca interacción entre los diferentes procesadores que forman el sistema. Atendiendo a esta y a otras características, la clasificación moderna divide a los sistemas en dos tipos: Sistemas multiprocesador (fuertemente acoplados) y sistemas multicomputadoras (débilmente acoplados). Multiprocesadores o Equipo Paralelo de Memoria Compartida Un multiprocesador puede verse como una computadora paralela compuesta por varios procesadores interconectados que comparten un mismo sistema de memoria. Los sistemas multiprocesadores son arquitecturas MIMD con memoria compartida. Tienen un único espacio de direcciones para todos los procesadores y los mecanismos de comunicación se basan en el paso de mensajes desde el punto de vista del programador. Dado que los multiprocesadores comparten diferentes módulos de memoria, pudiendo acceder a un mismo módulo varios procesadores, a los multiprocesadores también se les llama sistemas de memoria compartida. Figura 7: Arquitectura de una computadora paralela con memoria compartida Para hacer uso de la memoria compartida por más de un procesador, se requiere hacer uso de técnicas de semáforos que mantienen la integridad de la memoria; esta arquitectura no puede crecer mucho en el número de procesadores interconectados por la saturación rápida del bus o del medio de interconexión. Dependiendo de la forma en que los procesadores comparten la memoria, se clasifican en sistemas multiprocesador UMA, NUMA, COMA y Pipeline. 89 Uniform Memory Access (UMA) Sistema multiprocesador con acceso uniforme a memoria. La memoria física es uniformemente compartida por todos los procesadores, esto quiere decir que todos los procesadores tienen el mismo tiempo de acceso a todas las palabras de la memoria. Cada procesador tiene su propia caché privada y también se comparten los periféricos. Los multiprocesadores son sistemas fuertemente acoplados (tightly-coupled), dado el alto grado de compartición de los recursos (hardware o software) y el alto nivel de interacción entre procesadores, lo que hace que un procesador dependa de lo que hace otro. El sistema de interconexión debe ser rápido y puede ser de uno de los siguientes tipos: bus común, red crossbar y red multietapa. Este modelo es conveniente para aplicaciones de propósito general y de tiempo compartido por varios usuarios, existen dos categorías de sistemas UMA. Sistema Simétrico Cuando todos los procesadores tienen el mismo tiempo de acceso a todos los componentes del sistema (incluidos los periféricos), reciben el nombre de sistemas multiprocesador simétrico. Los procesadores tienen el mismo dominio (prioridad) sobre los periféricos y cada procesador tiene la misma capacidad para procesar. Sistema Asimétrico Los sistemas multiprocesador asimétrico, son sistemas con procesadores maestros y procesadores esclavos, en donde sólo los primeros pueden ejecutar aplicaciones y dónde en tiempo de acceso para diferentes procesadores no es el mismo. Los procesadores esclavos (attached) ejecutan código usuario bajo la supervisión del maestro, por lo tanto cuando una aplicación es ejecutada en un procesador maestro dispondrá de una cierta prioridad. Non Uniform Memory Access (NUMA) Un sistema multiprocesador NUMA es un sistema de memoria compartida donde el tiempo de acceso varía según donde se encuentre localizado el acceso. El acceso a memoria, por tanto, no es uniforme para diferentes procesadores, existen memorias locales asociadas a cada procesador y estos pueden acceder a datos de su memoria local de una manera más rápida que a las memorias de otros procesadores, debido a que primero debe aceptarse dicho acceso por el procesador del que depende el módulo de memoria local. Todas las memorias locales conforman la memoria global compartida y físicamente distribuida y accesible por todos los procesadores. Cache Only Memory Access (COMA) Los sistemas COMA son un caso especial de los sistemas NUMA. Este tipo de sistemas no ha tenido mucha trascendencia, al igual que los sistemas SIMD. 90 Las memorias distribuidas son memorias cachés, por este motivo es un sistema muy restringido en cuanto a la capacidad de memoria global. No hay jerarquía de memoria en cada módulo procesador. Todas las cachés forman un mismo espacio global de direcciones. El acceso a las cachés remotas se realiza a través de los directorios distribuidos de las cachés. Dependiendo de la red de interconexión utilizada, se pueden utilizar jerarquías en los directorios para ayudar a la localización de copias de bloques de caché. Procesador Vectorial Pipeline En la actualidad es común encontrar en un solo procesador los denominados Pipeline o Procesador Vectorial Pipeline del tipo MISD. En estos procesadores los vectores fluyen a través de las unidades aritméticas Pipeline. Las unidades constan de una cascada de etapas de procesamiento compuestas de circuitos que efectúan operaciones aritméticas o lógicas sobre el flujo de datos que pasan a través de ellas, las etapas están separadas por registros de alta velocidad usados para guardar resultados intermedios. Así la información que fluye entre las etapas adyacentes está bajo el control de un reloj que se aplica a todos los registros simultáneamente. Multicomputadoras o Equipo Paralelo de Memoria Distribuida Los sistemas multicomputadoras se pueden ver como una computadora paralela en el cual cada procesador tiene su propia memoria local. En estos sistemas la memoria se encuentra distribuida y no compartida como en los sistemas multiprocesador. Los procesadores se comunican a través de paso de mensajes, ya que éstos sólo tienen acceso directo a su memoria local y no a las memorias del resto de los procesadores. Figura 8: Arquitectura de una computadora paralela con memoria distribuida La transferencia de los datos se realiza a través de la red de interconexión que conecta un subconjunto de procesadores con otro subconjunto. La transferencia de unos procesadores a otros se realiza por múltiples transferencias entre procesadores conectados dependiendo del establecimiento de dicha red. 91 Dado que la memoria está distribuida entre los diferentes elementos de proceso, estos sistemas reciben el nombre de distribuidos. Por otra parte, estos sistemas son débilmente acoplados, ya que los módulos funcionan de forma casi independiente unos de otros. Este tipo de memoria distribuida es de acceso lento por ser peticiones a través de la red, pero es una forma muy efectiva de tener acceso a un gran volumen de memoria. Equipo Paralelo de Memoria Compartida-Distribuida La tendencia actual en las máquinas paralelas es de aprovechar las facilidades de programación que ofrecen los ambientes de memoria compartida y la escalabilidad de las ambientes de memoria distribuida. En este modelo se conectan entre si módulos de multiprocesadores, pero se mantiene la visión global de la memoria a pesar de que es distribuida. Clusters El desarrollo de sistemas operativos y compiladores del dominio público (Linux y software GNU), estándares para el pase de mensajes (MPI), conexión universal a periféricos (PCI), etc. han hecho posible tomar ventaja de los económicos recursos computacionales de producción masiva (CPU, discos, redes). La principal desventaja que presenta a los proveedores de multicomputadoras es que deben satisfacer una amplia gama de usuarios, es decir, deben ser generales. Esto aumenta los costos de diseños y producción de equipos, así como los costos de desarrollo de software que va con ellos: sistema operativo, compiladores y aplicaciones. Todos estos costos deben ser añadidos cuando se hace una venta. Por supuesto alguien que sólo necesita procesadores y un mecanismo de pase de mensajes no debería pagar por todos estos añadidos que nunca usará. Estos usuarios son los que están impulsando el uso de clusters principalmente de computadoras personales (PC), cuya arquitectura se muestra a continuación: Figura 9: Arquitectura de un cluster Los cluster se pueden clasificar en dos tipos según sus características físicas: Cluster homogéneo si todos los procesadores y/o nodos participantes en el equipo paralelo son iguales en capacidad de cómputo (en la cual es permitido variar la cantidad de memoria o disco duro en cada procesador). 92 Cluster heterogéneo es aquel en que al menos uno de los procesadores y/o nodos participantes en el equipo paralelo son de distinta capacidad de cómputo. Los cluster pueden formarse de diversos equipos; los más comunes son los de computadoras personales, pero es creciente el uso de computadoras multiprocesador de más de un procesador de memoria compartida interconectados por red con los demás nodos del mismo tipo, incluso el uso de computadoras multiprocesador de procesadores vectoriales Pipeline. Los cluster armados con la configuración anterior tienen grandes ventajas para procesamiento paralelo: • La reciente explosión en redes implica que la mayoría de los componentes necesarios para construir un cluster son vendidos en altos volúmenes y por lo tanto son económicos. Ahorros adicionales se pueden obtener debido a que sólo se necesitará una tarjeta de vídeo, un monitor y un teclado por cluster. El mercado de los multiprocesadores es más reducido y más costoso. • Remplazar un componente defectuoso en un cluster es relativamente trivial comparado con hacerlo en un multiprocesador, permitiendo una mayor disponibilidad de clusters cuidadosamente diseñados. Desventajas del uso de clusters de computadoras personales para procesamiento paralelo: • Con raras excepciones, los equipos de redes generales producidos masivamente no están diseñados para procesamiento paralelo y típicamente su latencia es alta y los anchos de banda pequeños comparados con multiprocesadores. Dado que los clusters explotan tecnología que sea económica, los enlaces en el sistema no son veloces implicando que la comunicación entre componentes debe pasar por un proceso de protocolos de negociación lentos, incrementando seriamente la latencia. En muchos y en el mejor de los casos (debido a costos) se recurre a una red tipo Fast Ethernet restringimiento la escalabilidad del cluster. • Hay poco soporte de software para manejar un cluster como un sistema integrado. • Los procesadores no son tan eficientes como los procesadores usados en los multiprocesadores para manejar múltiples usuarios y/o procesos. Esto hace que el rendimiento de los clusters se degrade con relativamente pocos usuarios y/o procesos. • Muchas aplicaciones importantes disponibles en multiprocesadores y optimizadas para ciertas arquitecturas, no lo están en clusters. Sin lugar a duda los clusters presentan una alternativa importante para varios problemas particulares, no sólo por su economía, si no también porque 93 pueden ser diseñados y ajustados para ciertas aplicaciones. Las aplicaciones que pueden sacar provecho de clusters son en donde el grado de comunicación entre procesos es de bajo a medio. Tipos de Cluster Básicamente existen tres tipo de clusters, cada uno de ellos ofrece ventajas y desventajas, el tipo más adecuado para el cómputo científico es del de altorendimiento, pero existen aplicaciones científicas que pueden usar más de un tipo al mismo tiempo. • Alta-disponibilidad (Fail-over o High-Availability): este tipo de cluster esta diseñado paramantener uno o varios servicios disponibles incluso a costa de rendimiento, ya que su función principal es que el servicio jamás tenga interrupciones como por ejemplo un servicio de bases de datos. • Alto-rendimiento (HPC o High Performance Computing): este tipo de cluster está diseñado para obtener el máximo rendimiento de la aplicación utilizada incluso a costa de la disponibilidad del sistema, es decir el cluster puede sufrir caídas, este tipo de configuración está orientada a procesos que requieran mucha capacidad de cálculo. • Balanceo de Carga (Load-balancing): este tipo de cluster esta diseñado para balancear la carga de trabajo entre varios servidores, lo que permite tener, por ejemplo, un servicio de cálculo intensivo multiusuarios que detecte tiempos muertos del proceso de un usuario para ejecutar en dichos tiempos procesos de otros usuarios. Grids Son cúmulos (grupo de clusters) de arquitecturas en paralelo interconectados por red, los cuales distribuyen tareas entre los clusters que lo forman, estos pueden ser homogéneos o heterogéneos en cuanto a los nodos componentes del cúmulo. Este tipo de arquitecturas trata de distribuir cargas de trabajo acorde a las características internas de cada cluster y las necesidades propias de cada problema, esto se hace a dos niveles, una en la parte de programación conjuntamente con el balance de cargas y otra en la parte de hardware que tiene que ver con las características de cada arquitectura que conforman al cúmulo. 7.2. Métricas de Desempeño Las métricas de desempeño del procesamiento de alguna tarea en paralelo es un factor importante para medir la eficiencia y consumo de recursos al resolver una tarea con un número determinado de procesadores y recursos relacionados de la interconexión de éstos. Entre las métricas para medir desempeño en las cuales como premisa se mantiene fijo el tamaño del problema, destacan las siguientes: Factor de aceleración, eficiencia y fracción serial. Cada una de ellas mide algo en particular 94 y sólo la combinación de estas dan un panorama general del desempeño del procesamiento en paralelo de un problema en particular en una arquitectura determinada al ser comparada con otras. Factor de Aceleración (o Speed-Up) Se define como el cociente del tiempo que se tarda en completar el cómputo de la tarea usando un sólo procesador entre el tiempo que necesita para realizarlo en p procesadores trabajando en paralelo s = T (1) T (p) (314) en ambos casos se asume que se usará el mejor algoritmo tanto para un solo procesador como para p procesadores. Esta métrica en el caso ideal debería de aumentar de forma lineal al aumento del número de procesadores. Eficiencia Se define como el cociente del tiempo que se tarda en completar el cómputo de la tarea usando un solo procesador entre el número de procesadores multiplicado por el tiempo que necesita para realizarlo en p procesadores trabajando en paralelo e = T (1) pT (p) = s p . (315) Este valor será cercano a la unidad cuando el hardware se esté usando de manera eficiente, en caso contrario el hardware será desaprovechado. Fracción serial Se define como el cociente del tiempo que se tarda en completar el cómputo de la parte secuencial de una tarea entre el tiempo que se tarda el completar el cómputo de la tarea usando un solo procesador f = Ts T (1) (316) pero usando la ley de Amdahl T (p) = Ts + Tp p y rescribiéndola en términos de factor de aceleración, obtenemos la forma operativa del cálculo de la fracción serial que adquiere la forma siguiente f = 1 s − 1 p 1 − 1 p . (317) Esta métrica permite ver las inconsistencias en el balance de cargas, ya que su valor debiera de tender a cero en el caso ideal, por ello un incremento en el valor de f es un aviso de granularidad fina con la correspondiente sobrecarga en los procesos de comunicación. 95 7.3. Cómputo Paralelo para Sistemas Continuos Como se mostró en los capítulos anteriores, la solución de los sistemas continuos usando ecuaciones diferenciales parciales genera un alto consumo de memoria e involucran un amplio tiempo de procesamiento; por ello nos interesa trabajar en computadoras que nos puedan satisfacer estas demandas. Actualmente, en muchos centros de cómputo es una práctica común usar directivas de compilación en equipos paralelos sobre programas escritos de forma secuencial, con la esperanza que sean puestos por el compilador como programas paralelos. Esto en la gran mayoría de los casos genera códigos poco eficientes, pese a que corren en equipos paralelos y pueden usar toda la memoria compartida de dichos equipos, el algoritmo ejecutado continua siendo secuencial en la gran mayoría del código. Si la arquitectura paralela donde se implemente el programa es UMA de acceso simétrico, los datos serán accesados a una velocidad de memoria constante. En caso contrario, al acceder a un conjunto de datos es común que una parte de estos sean locales a un procesador (con un acceso del orden de nano segundos), pero el resto de los datos deberán de ser accesados mediante red (con acceso del orden de mili segundos), siendo esto muy costoso en tiempo de procesamiento. Por ello, si usamos métodos de descomposición de dominio es posible hacer que el sistema algebraico asociado pueda distribuirse en la memoria local de múltiples computadoras y que para encontrar la solución al problema se requiera poca comunicación entre los procesadores. Por lo anterior, si se cuenta con computadoras con memoria compartida o que tengan interconexión por bus, salvo en casos particulares no será posible explotar éstas características eficientemente. Pero en la medida en que se adecuen los programas para usar bibliotecas y compiladores acordes a las características del equipo disponible (algunos de ellos sólo existen de manera comercial) la eficiencia aumentará de manera importante. La alternativa más adecuada (en costo y flexibilidad), es trabajar con computadoras de escritorio interconectadas por red que pueden usarse de manera cooperativa para resolver nuestro problema. Los gastos en la interconexión de los equipos son mínimos (sólo el switch y una tarjeta de red por equipo y cables para su conexión). Por ello los clusters y los grids son en principio una buena opción para la resolución de este tipo de problemas. Esquema de Paralelización Maestro-Esclavo La implementación de los métodos de descomposición de dominio que se trabajarán será mediante el esquema Maestro-Esclavo (Farmer) en el lenguaje de programación C++ bajo la interfaz de paso de mensajes MPI trabajando en un cluster Linux Debian. Donde tomando en cuenta la implementación en estrella del cluster, el modelo de paralelismo de MPI y las necesidades propias de comunicación del programa, el nodo maestro tendrá comunicación sólo con cada nodo esclavo y no existirá comunicación entre los nodos esclavos, esto reducirá las comunicaciones y optimizará el paso de mensajes. El esquema de paralelización maestro-esclavo, permite sincronizar por parte 96 del nodo maestro las tareas que se realizan en paralelo usando varios nodos esclavos, éste modelo puede ser explotado de manera eficiente si existe poca comunicación entre el maestro y el esclavo y los tiempos consumidos en realizar las tareas asignadas son mayores que los períodos involucrados en las comunicaciones para la asignación de dichas tareas. De esta manera se garantiza que la mayoría de los procesadores estarán trabajando de manera continua y existirán pocos tiempos muertos. Figura 10: Esquema del maestro-esclavo Un factor limitante en este esquema es que el nodo maestro deberá de atender todas las peticiones hechas por cada uno de los nodos esclavos, esto toma especial relevancia cuando todos o casi todos los nodos esclavos compiten por ser atendidos por el nodo maestro. Se recomienda implementar este esquema en un cluster heterogéneo en donde el nodo maestro sea más poderoso computacionalmente que los nodos esclavos. Si a éste esquema se le agrega una red de alta velocidad y de baja latencia, se le permitirá operar al cluster en las mejores condiciones posibles, pero este esquema se verá degradado al aumentar el número de nodos esclavos inexorablemente. Pero hay que ser cuidadosos en cuanto al número de nodos esclavos que se usan en la implementación en tiempo de ejecución versus el rendimiento general del sistema al aumentar estos, algunas observaciones posibles son: • El esquema maestro-esclavo programado en C++ y usando MPI lanza P procesos (uno para el nodo maestro y P − 1 para los nodos esclavos), estos en principio corren en un solo procesador pero pueden ser lanzados en múltiples procesadores usando una directiva de ejecución, de esta manera es posible que en una sola maquina se programe, depure y sea puesto a punto el código usando mallas pequeñas (del orden de cientos de nodos) y cuando este listo puede mandarse a producción en un cluster. • El esquema maestro-esclavo no es eficiente si sólo se usan dos procesadores (uno para el nodo maestro y otro para el nodo esclavo), ya que el nodo maestro en general no realiza los cálculos pesados y su principal función será la de distribuir tareas; los cálculos serán delegados al nodo esclavo. En el caso que nos interesa implementar, el método de descomposición de dominio adolece de este problema. 97 Paso de Mensajes Usando MPI Para poder intercomunicar al nodo maestro con cada uno de los nodos esclavos se usa la interfaz de paso de mensajes (MPI), una biblioteca de comunicación para procesamiento en paralelo. MPI ha sido desarrollado como un estándar para el paso de mensajes y operaciones relacionadas. Este enfoque es adoptado por usuarios e implementadores de bibliotecas, en la cual se proveen a los programas de procesamiento en paralelo de portabilidad y herramientas necesarias para desarrollar aplicaciones que puedan usar el cómputo paralelo de alto desempeño. El modelo de paso de mensajes posibilita a un conjunto de procesos que tienen solo memoria local la comunicación con otros procesos (usando Bus o red) mediante el envío y recepción de mensajes. Por definición el paso de mensajes posibilita transferir datos de la memoria local de un proceso a la memoria local de cualquier otro proceso que lo requiera. En el modelo de paso de mensajes para equipos paralelos, los procesos se ejecutan en paralelo, teniendo direcciones de memoria separada para cada proceso, la comunicación ocurre cuando una porción de la dirección de memoria de un proceso es copiada mediante el envío de un mensaje dentro de otro proceso en la memoria local mediante la recepción del mismo. Las operaciones de envío y recepción de mensajes es cooperativa y ocurre sólo cuando el primer proceso ejecuta una operación de envío y el segundo proceso ejecuta una operación de recepción, los argumentos base de estas funciones son: • Para el que envía, la dirección de los datos a transmitir y el proceso destino al cual los datos se enviarán. • Para el que recibe, debe de tener la dirección de memoria donde se pondrán los datos recibidos, junto con la dirección del proceso del que los envío. Es decir: Send(dir, lg, td, dest, etiq, com) {dir, lg, td} describe cuántas ocurrencias lg de elementos del tipo de dato td se transmitirán empezando en la dirección de memoria dir. {des, etiq, com} describe el identificador etq de destino des asociado con la comunicación com. Recv(dir, mlg, td, fuent, etiq, com, st) {dir, lg, td} describe cuántas ocurrencias lg de elementos del tipo de dato td se transmitirán empezando en la dirección de memoria dir. {f uent, etiq, com, est} describe el identificador etq de la fuente fuent asociado con la comunicación com y el estado st. El conjunto básico de directivas (en nuestro caso sólo se usan estas) en C++ de MPI son: 98 MPI::Init Inicializa al MPI MPI::COMM_WORLD.Get_size Busca el número de procesos existentes MPI::COMM_WORLD.Get_rank Busca el identificador del proceso MPI::COMM_WORLD.Send Envía un mensaje MPI::COMM_WORLD.Recv Recibe un mensaje MPI::Finalize Termina al MPI Estructura del Programa Maestro-Esclavo La estructura del programa se realizo para que el nodo maestro mande trabajos de manera síncrona a los nodos esclavos. Cuando los nodos esclavos terminan la tarea asignada, avisan al nodo maestro para que se le asigne otra tarea (estas tareas son acordes a la etapa correspondiente del método de descomposición de dominio ejecutándose en un instante dado). En la medida de lo posible se trata de mandar paquetes de datos a cada nodo esclavo y que estos regresen también paquetes al nodo maestro, a manera de reducir las comunicaciones al mínimo y tratar de mantener siempre ocupados a los nodos esclavos para evitar los tiempos muertos, logrando con ello una granularidad gruesa, ideal para trabajar con clusters. La estructura básica del programa bajo el esquema maestro-esclavo codificada en C++ y usando MPI es: main(int argc, char *argv[]) { MPI::Init(argc,argv); ME_id = MPI::COMM_WORLD.Get_rank(); MP_np = MPI::COMM_WORLD.Get_size(); if (ME_id == 0) { // Operaciones del Maestro } else { // Operaciones del esclavo con identificador ME_id } MPI::Finalize(); } En este único programa se deberá de codificar todas las tareas necesarias para el nodo maestro y cada uno de los nodos esclavos, así como las formas de intercomunicación entre ellos usando como distintivo de los distintos procesos a la variable ME_id. Para más detalles de esta forma de programación y otras funciones de MPI ver [21] y [6]. Los factores limitantes para el esquema maestro-esclavo pueden ser de dos tipos, los inherentes al propio esquema maestro-esclavo y al método de descomposición de dominio: 99 • El esquema de paralelización maestro-esclavo presupone contar con un nodo maestro lo suficientemente poderoso para atender simultáneamente las tareas síncronas del método de descomposición de dominio, ya que este distribuye tareas acorde al número de subdominios, estas si son balanceadas ocasionaran que todos los procesadores esclavos terminen al mismo tiempo y el nodo maestro tendrá que atender múltiples comunicaciones simultáneamente, degradando su rendimiento al aumentar el número de subdominios. • Al ser síncrono el método de descomposición de dominio, si un nodo esclavo acaba la tarea asignada y avisa al nodo maestro, este no podrá asignarle otra tarea hasta que todos los nodos esclavos concluyan la suya. Para los factores limitantes inherente al propio esquema maestro-esclavo, es posible implementar algunas operaciones del nodo maestro en paralelo, ya sea usando equipos multiprocesador o en más de un nodo distintos a los nodos esclavos. Para la parte inherente al método de descomposición de dominio, la parte medular la da el balanceo de cargas. Es decir que cada nodo esclavo tenga una carga de trabajo igual al resto de los nodos. Este balanceo de cargas puede no ser homogéneo por dos razones: Al tener P procesadores en el equipo paralelo, la descomposición del dominio no sea la adecuada. Si se tiene una descomposición particular, esta se implemente en un número de procesadores inadecuado. Cualquiera de las dos razones generarán desbalanceo de la carga en los nodos esclavos, ocasionando una perdida de eficiencia en el procesamiento de un problema bajo una descomposición particular en una configuración del equipo paralelo especifica, es por esto que en algunos casos al aumentar el número de procesadores que resuelvan la tarea no se aprecia una disminución del de tiempo de procesamiento. El número de procesadores P que se usen para resolver un dominio Ω y tener buen balance de cargas puede ser conocido si aplicamos el siguiente procedimiento: Si el dominio Ω se descompone en n ×m subdominios (la partición gruesa), entonces se generarán s = n ∗ m subdominios Ωi, en este caso, se tiene un buen balanceo de cargas si (P − 1) | s. La partición fina se obtiene al descomponer a cada subdominio Ωi en p × q subdominios. Como ejemplo, supongamos que deseamos resolver el dominio Ω usando 81× 81 nodos (nodos = n∗p+1 y nodos = m∗q+1), de manera inmediata nos surgen las siguientes preguntas: ¿cuales son las posibles descomposiciones validas? y ¿en cuantos procesadores se pueden resolver cada descomposición?. Para este ejemplo, sin hacer la tabla exhaustiva obtenemos: 100 Partición Subdominios Procesadores 1x2 y 80x40 2 2,3 1x4 y 80x20 5 2,5 1x5 y 80x16 6 2,6 2x1 y 40x80 2 2,3 2x2 y 40x40 4 2,3,5 2x4 y 40x20 8 2,3,5,9 2x5 y 40x16 10 2,3,6,11 2x8 y 40x10 16 2,3,5,9,17 4x1 y 20x80 4 2,3,5 4x2 y 20x40 8 2,3,5,9 4x4 y 20x20 16 2,3,5,9,17 4x5 y 20x16 20 2,3,5,6,11,21 5x1 y 16x80 5 2,6 5x2 y 16x40 10 2,3,6,11 5x4 y 16x20 20 2,3,5,6,11,21 5x5 y 16x16 25 2,6,26 De esta tabla es posible seleccionar (para este ejemplo en particular), las descomposiciones que se adecuen a las necesidades particulares del equipo con que se cuente. Sin embargo hay que tomar en cuenta siempre el número de nodos por subdominio de la partición fina, ya que un número de nodos muy grande puede que exceda la cantidad de memoria que tiene el nodo esclavo y un número pequeño estaría infrautilizando el poder computacional de los nodos esclavos. De las particiones seleccionadas se pueden hacer corridas de prueba para evaluar su rendimiento, hasta encontrar la que menor tiempo de ejecución consuma, maximizando así la eficiencia del equipo paralelo. Programación Paralela en Multihilos En una computadora, sea secuencial o paralela, para aprovechar las capacidades crecientes del procesador, el sistema operativo divide su tiempo de procesamiento entre los distintos procesos, de forma tal que para poder ejecutar a un proceso, el kernel les asigna a cada proceso una prioridad y con ello una fracción del tiempo total de procesamiento, de forma tal que se pueda atender a todos y cada uno de los procesos de manera eficiente. En particular, en la programación en paralelo usando MPI, cada proceso (que eventualmente puede estar en distinto procesador) se lanza como una copia del programa con datos privados y un identificador del proceso único, de tal forma que cada proceso sólo puede compartir datos con otro proceso mediante paso de mensajes. Esta forma de lanzar procesos por cada tarea que se desee hacer en paralelo es costosa, por llevar cada una de ellas todo una gama de subprocesos para poderle asignar recursos por parte del sistema operativo. Una forma más eficiente de hacerlo es que un proceso pueda generar bloques de subprocesos que puedan ser ejecutados como parte del proceso (como subtareas), así en el tiempo asignado 101 se pueden atender a más de un subproceso de manera más eficiente, esto es conocido como programación multihilos. Los hilos realizarán las distintas tareas necesarias en un proceso. Para hacer que los procesos funcionen de esta manera, se utilizan distintas técnicas que le indican kernel cuales son las partes del proceso que pueden ejecutarse simultáneamente y el procesador asignará una fracción de tiempo exclusivo al hilo del tiempo total asignado al proceso. Los datos pertenecientes al proceso pasan a ser compartidos por los subprocesos lanzados en cada hilo y mediante una técnica de semáforos el kernel mantiene la integridad de estos. Esta técnica de programación puede ser muy eficiente si no se abusa de este recurso, permitiendo un nivel más de paralelización en cada procesador. Esta forma de paralelización no es exclusiva de equipos multiprocesadores o multicomputadoras, ya que pueden ser implementados a nivel de sistema operativo. 102 8. Implementación Computacional Secuencial y Paralela de DDM A partir de los modelos matemáticos (capítulo 2 y 3) y los modelos numéricos (capítulos 4, 5, y 6), en este capítulo se describe el modelo computacional contenido en un programa de cómputo orientado a objetos en el lenguaje de programación C++ en su forma secuencial y en su forma paralela en C++ usando la interfaz de paso de mensajes (MPI) bajo el esquema maestro-esclavo (capítulo 7). Esto no sólo nos ayudará a demostrar que es factible la construcción del propio modelo computacional a partir del modelo matemático y numérico para la solución de problemas reales. Además, se mostrará los alcances y limitaciones en el consumo de los recursos computacionales, evaluando algunas de las variantes de los métodos numéricos con los que es posible implementar el modelo computacional y haremos el análisis de rendimiento sin llegar a ser exhaustivo esté. También exploraremos los alcances y limitaciones de cada uno de los métodos implementados (FEM, DDM secuencial y paralelo) y como es posible optimizar los recursos computacionales con los que se cuente. Primeramente hay que destacar que el paradigma de programación orientada a objetos es un método de implementación de programas, organizados como colecciones cooperativas de objetos. Cada objeto representa una instancia de alguna clase y cada clase es miembro de una jerarquía de clases unidas mediante relaciones de herencia, contención, agregación o uso. Esto nos permite dividir en niveles la semántica de los sistemas complejos tratando así con las partes, que son más manejables que el todo, permitiendo su extensión y un mantenimiento más sencillo. Así, mediante la herencia, contención, agregación o usó nos permite generar clases especializadas que manejan eficientemente la complejidad del problema. La programación orientada a objetos organiza un programa entorno a sus datos (atributos) y a un conjunto de interfases bien definidas para manipular estos datos (métodos dentro de clases reusables) esto en oposición a los demás paradigmas de programación. El paradigma de programación orientada a objetos sin embargo sacrifica algo de eficiencia computacional por requerir mayor manejo de recursos computacionales al momento de la ejecución. Pero en contraste, permite mayor flexibilidad al adaptar los códigos a nuevas especificaciones. Adicionalmente, disminuye notoriamente el tiempo invertido en el mantenimiento y búsqueda de errores dentro del código. Esto tiene especial interés cuando se piensa en la cantidad de meses invertidos en la programación comparado con los segundos consumidos en la ejecución del mismo. Para empezar con la implementación computacional, primeramente definiremos el problema a trabajar. Este, pese a su sencillez, no pierde generalidad permitiendo que el modelo mostrado sea usado en muchos sistemas de la ingeniería y la ciencia. 103 8.1. El Operador de Laplace y la Ecuación de Poisson Consideramos como modelo matemático el problema de valor en la frontera (BVP) asociado con el operador de Laplace en dos dimensiones, el cual en general es usualmente referido como la ecuación de Poisson, con condiciones de frontera Dirichlet, definido en Ω como: −∇2u = fΩ en Ω (318) u = g∂Ω en ∂Ω. Se toma está ecuación para facilitar la comprensión de las ideas básicas. Es un ejemplo muy sencillo, pero gobierna los modelos de muchos sistemas de la ingeniería y de la ciencia, entre ellos el flujo de agua subterránea a través de un acuífero isotrópico, homogéneo bajo condiciones de equilibrio y es muy usada en múltiples ramas de la física. Por ejemplo, gobierna la ecuación de la conducción de calor en un sólido bajo condiciones de equilibrio. En particular consideramos el problema con Ω definido en: Ω = [−1, 1] × [0, 1] (319) donde fΩ = 2n2π2 sin(nπx) ∗ sin(nπy) y g∂Ω = 0 (320) cuya solución es u(x, y) = sin(nπx) ∗ sin(nπy). (321) Para las pruebas de rendimiento en las cuales se evalúa el desempeño de los programas realizados se usa n = 10, pero es posible hacerlo con n ∈ N grande. Por ejemplo para n = 4, la solución es u(x, y) = sin(4πx)∗sin(4πy), cuya gráfica se muestra a continuación: Figura 11: Solución a la ecuación de Poisson para n=4. 104 Hay que hacer notar que al implementar la solución numérica por el método del elemento finito y el método de subestructuración secuencial en un procesador, un factor limitante para su operación es la cantidad de memoria disponible en la computadora, ya que el sistema algebraico de ecuaciones asociado a esté problema crece muy rápido (del orden de n2), donde n es el número de nodos en la partición. Es por ello que la elección de un buen manejador de matrices será determinante en la eficiencia alcanzada por las distintas implementaciones de los programas. Actualmente existen múltiples bibliotecas que permiten manipular operaciones de matrices tanto en forma secuencial como en paralelo (hilos y Pipeline) para implementarlas tanto en procesadores con memoria compartida como distribuida. Pero no están presentes en todas las arquitecturas y sistemas operativos. Por ello en este trabajo se implementaron todas las operaciones necesarias usando clases que fueron desarrolladas sin usar ninguna biblioteca externa a las proporcionadas por el compilador de C++ de GNU, permitiendo la operación de los programas desarrollados en múltiples sistemas operativos del tipo Linux, Unix y Windows. En cuanto a las pruebas de rendimiento que mostraremos en las siguientes secciones se usaron para la parte secuencial el equipo: • Computadora Pentium IV HT a 2.8 GHtz con 1 GB de RAM corriendo bajo el sistema operativo Linux Debian Stable con el compilador g++ de GNU. Para la parte paralela se usaron los equipos siguientes: • Cluster homogéneo de 10 nodos duales Xeon a 2.8 GHtz con 1 GB de RAM por nodo, unidos mediante una red Ethernet de 1 Gb, corriendo bajo el sistema operativo Linux Debian Stable con el compilador mpiCC de MPI de GNU. • Cluster heterogéneo con el nodo maestro Pentium IV HT a 3.4 GHtz con 1 GB de RAM y 7 nodos esclavos Pentium IV HT a 2.8 GHtz con 0.5 GB de RAM por nodo, unidos mediante una red Ethernet de 100 Mb, corriendo bajo el sistema operativo Linux Debian Stable con el compilador mpiCC de MPI de GNU. A estos equipos nos referiremos en lo sucesivo como equipo secuencial, cluster homogéneo y cluster heterogéneo respectivamente. El tiempo dado en los resultados de las distintas pruebas de rendimiento de los programas y mostrado en todas las tablas y gráficas fue tomado como un promedio entre por lo menos 5 corridas, redondeado el resultado a la unidad siguiente. En todos los cálculos de los métodos numéricos usados para resolver el sistema lineal algebraico asociado se usó una tolerancia mínima de 1×10−10. Ahora, veremos la implementación del método de elemento finito secuencial para después continuar con el método de descomposición de dominio tanto secuencial como paralelo y poder analizar en cada caso los requerimientos de cómputo, necesarios para correr eficientemente un problema en particular. 105 8.2. Método del Elemento Finito Secuencial A partir de la formulación del método de elemento finito visto en la sección (4.2), la implementación computacional que se desarrolló tiene la jerarquía de clases siguiente: Figura 12: Jerarquía de clases para el método de elemento finito Donde las clases participantes en FEM2D Rectángulos son: La clase Interpolador Lineal define los interpoladores lineales usados por el método de elemento finito. La clase Problema define el problema a tratar, es decir, la ecuación diferencial parcial, valores de frontera y dominio. La clase Base FEM ayuda a definir los nodos al usar la clase Geometría y mantiene las matrices generadas por el método y a partir de la clase Resuelve Ax=B se dispone de diversas formas de resolver el sistema lineal asociado al método. La clase FEM2D controla lo necesario para poder hacer uso de la geometría en 2D y conocer los nodos interiores y de frontera, con ellos poder montar la matriz de rigidez y ensamblar la solución. La clase FEM2D Rectángulos permite calcular la matriz de rigidez para generar el sistema algebraico de ecuaciones asociado al método. Notemos que esta misma jerarquía permite trabajar problemas en una y dos dimensiones, en el caso de dos dimensiones podemos discretizar usando 106 rectángulos o triángulos, así como usar varias opciones para resolver el sistema lineal algebraico asociado a la solución de EDP. Como ya se menciono, el método de elemento finito es un algoritmo secuencial, por ello se implementa para que use un solo procesador y un factor limitante para su operación es la cantidad de memoria disponible en la computadora, por ejemplo: Resolver la Ec. (318) con una partición rectangular de 513×513 nodos, genera 262144 elementos rectangulares con 263169 nodos en total, donde 261121 son desconocidos; así el sistema algebraico de ecuaciones asociado a este problema es de dimensión 261121 × 261121. Usando el equipo secuencial, primeramente evaluaremos el desempeño del método de elemento finito con los distintos métodos para resolver el sistema algebraico de ecuaciones, encontrando los siguientes resultados: Método Iterativo Iteraciones Tiempo Total Jacobi 865037 115897 seg. Gauss-Seidel 446932 63311 seg. Gradiente Conjugado 761 6388 seg. Como se observa el uso del método de gradiente conjugado es por mucho la mejor elección. En principio, podríamos quedarnos solamente con el método de gradiente conjugado sin hacer uso de precondicionadores por los buenos rendimientos encontrados hasta aquí, pero si se desea resolver un problema con un gran número de nodos, es conocido el aumento de eficiencia al hacer uso de precondicionadores. Ahora, si tomamos ingenuamente el método de elemento finito conjuntamente con el método de gradiente conjugado con precondicionadores a posteriori (los más sencillos de construir) para resolver el sistema algebraico de ecuaciones, encontraremos los siguientes resultados: Precondicionador Iteraciones Tiempo Total Jacobi 760 6388 seg. SSOR 758 6375 seg. Factorización Incompleta 745 6373 seg. Como es notorio el uso del método de gradiente conjugado precondicionado con precondicionadores a posteriori no ofrece una ventaja significativa que compense el esfuerzo computacional invertido al crear y usar un precondicionador en los cálculos por el mal condicionamiento del sistema algebraico. Existen también precondicionadores a priori para el método de elemento finito, pero no es costeable en rendimiento su implementación. Finalmente, para el método de elemento finito las posibles mejoras de eficiencia para disminuir el tiempo de ejecución pueden ser: • Al momento de compilar los códigos usar directivas de optimización (ofrece mejoras de rendimiento en ejecución de 30% aproximadamente en las pruebas realizadas). 107 • Usar la biblioteca Lapack++ de licencia GNU que optimiza las operaciones en el manejo de los elementos de la matriz usando punteros y hacen uso de matrices bandadas (obteniéndose una mejora del rendimiento de 15% aproximadamente en las pruebas realizadas). • Combinando las opciones anteriores se obtiene una mejora sustancial de rendimiento en la ejecución (de 45% aproximadamente en las pruebas realizadas). Adicionalmente si se cuenta con un equipo con más de un procesador con memoria compartida es posible usar bibliotecas para la manipulación de matrices y vectores que paralelizan o usan Pipeline como una forma de mejorar el rendimiento del programa. Este tipo de mejoras cuando es posible usarlas disminuyen sustancialmente el tiempo de ejecución, ya que en gran medida el consumo total de CPU está en la manipulación de matrices, pero esto no hace paralelo al método de elemento finito. 8.3. Método de Subestructuración Secuencial A partir de la formulación del método de subestructuración visto en la sección podría ser construida sumando las S bi respectivamente según el orden de los nodos globales versus los nodos locales a cada subdominio. El sistema lineal virtual resultante SuΣ = b (323) es resuelto usando el método de gradiente conjugado visto en la sección (5.2), para ello no es necesario construir la matriz S con las contribuciones de cada Si correspondientes al subdominio i. Lo que hacemos es pasar a cada subdominio el vector uΣ i correspondiente a la i-ésima iteración del método de gradiente conjugado para que en cada subdominio se evalué ˜uΣ i = S i uΣ i localmente y con el resultado se forma el vector ˜uΣ = PE i=1 ˜uΣ i y se continué con los demás pasos del método. La implementación computacional que se desarrolló tiene una jerarquía de clases en la cual se agregan las clases FEM2D Rectángulos y Geometría, además de heredar a la clase Problema. De esta forma se rehusó todo el código desarrollado para FEM2D Rectángulos, la jerarquía queda como: 108 La clase DDM2D realiza la partición gruesa del dominio mediante la clase Geometría y controla la partición de cada subdominio mediante un objeto de la clase de FEM2D Rectángulos generando la partición fina del dominio. La resolución de los nodos de la frontera interior se hace mediante el método de gradiente conjugado, necesaria para resolver los nodos internos de cada subdominio. Figura 13: Jerarquía de clases para el método de subestructuración secuencial Así, el dominio Ω es descompuesto en una descomposición gruesa de n × m subdominios y cada subdominio Ωi se parte en p×q subdominios, generando la participación fina del dominio como se muestra en la figura: Figura 14: Descomposición del dominio Ω en E = n × m subdominios y cada subdominio Ωi en p × q subdominios El método de descomposición de dominio se implementó realizando las siguientes tareas: 109 A) La clase DDM2D genera la descomposición gruesa del dominio mediante la agregación de un objeto de la clase Geometría (supongamos particionado en n × m subdominios, generando s = n ∗ m subdominios Ωi, i = 1, 2, ..., E). B) Con esa geometría se construyen los objetos de FEM2D Rectángulos (uno por cada subdominio Ωi), donde cada subdominio es particionado (supongamos en p×q subdominios) y regresando las coordenadas de los nodos de frontera del subdominio correspondiente a la clase DDM2D. C) Con estas coordenadas, la clase DDM2D conoce a los nodos de la frontera interior (son estos los que resuelve el método de descomposición de dominio). Las coordenadas de los nodos de la frontera interior se dan a conocer a los objetos FEM2D Rectángulos, transmitiendo sólo aquellos que están en su subdominio. D) Después de conocer los nodos de la frontera interior, cada objeto FEM2D Rectángulos calcula las matrices AΣΣ sin realizar comunicación alguna. Al terminar de calcular las matrices se avisa a la clase DDM2D de la finalización de los cálculos. E) Mediante la comunicación de vectores del tamaño del número de nodos de la frontera interior entre la clase DDM2D y los objetos FEM2D Rectángulos, se prepara todo lo necesario para empezar el hmPétodo de gradiente conjugado y resolver el sistema lineal virtual E i=1 S i i uΣ = hPE i=1 bi i . F) Para usar el método de gradiente conjugado, se transmite un vector del tamaño del número de nodos de la frontera interior para que en cada objeto se realicen las operaciones pertinentes y resolver así el sistema algebraico asociado, esta comunicación se realiza de ida y vuelta entre la clase DDM2D y los objetos FEM2D Rectángulos tantas veces como iteraciones haga el método. Resolviendo con esto los nodos de la frontera interior uΣi. G) Al término de las iteraciones se pasa la solución uΣi de los nodos de la frontera interior que pertenecen a cada subdominio dentro de cada objeto FEM2D Rectángulos para que se resuelvan los nodos interiores uI i = ³ AII i ´ −1 ³ bI i − AIΣ i uΣi ´ , sin realizar comunicación alguna en el proceso, al concluir se avisa a la clase DDM2D de ello. I) La clase DDM2D mediante un último mensaje avisa que se concluya el programa, terminado así el esquema maestro-esclavo secuencial. 110 Por ejemplo, para resolver la Ec. (318), usando 513 × 513 nodos (igual al ejemplo de FEM2D Rectángulos secuencial), podemos tomar alguna de las siguientes descomposiciones: D es com p o s ición N odos In t e r io r e s S ub dominios E lem ento s Su b dom inio To tal N o do s S ub dominio N odos D escono cidos S ub dominio 2x2 y 256x256 260100 4 65536 66049 65025 4x4 y 128x128 258064 16 16384 16641 16129 8x8 y 64x64 254016 64 4096 4225 3969 16x16 y 32x32 246016 256 1024 1089 961 32x32 y 16x16 230400 1024 256 289 225 Cada una de las descomposiciones genera un problema distinto. Usando el equipo secuencial y evaluando el desempeño del método de subestructuración secuencial se obtuvieron los siguientes resultados: Partición Nodos Frontera Interior Iteraciones Tiempo Total 2x2 y 256x256 1021 139 5708 seg. 4x4 y 128x128 3057 159 2934 seg. 8x8 y 64x64 7105 204 1729 seg. 16x16 y 32x32 15105 264 1077 seg. 32x32 y 16x16 30721 325 1128 seg Nótese que aún en un solo procesador es posible encontrar una descomposición que desminuya los tiempos de ejecución (la descomposición de 2x2 y 256x256 concluye en 5708 seg. versus los 6388 seg. en el caso de FEM2D Rectángulos), ello es debido a que al descomponer el dominio en múltiples subdominios, la complejidad del problema es también disminuida y esto se ve reflejado en la disminución del tiempo de ejecución. En la última descomposición, en lugar de disminuir el tiempo de ejecución este aumenta, esto se debe a que se construyen muchos objetos FEM2D Rectángulos (1024 en este caso), con los cuales hay que hacer comunicación resultando muy costoso computacionalmente. Finalmente las posibles mejoras de eficiencia para el método de subestructuración secuencial para disminuir el tiempo de ejecución son las mismas que en el caso del método de elemento finito pero además se tienen que: • Encontrar la descomposición pertinente entre las posibles descomposiciones que consuma el menor tiempo de cálculo. Adicionalmente si se cuenta con un equipo con más de un procesador con memoria compartida es posible usar bibliotecas para la manipulación de matrices y vectores que paralelizan o usan Pipeline como una forma de mejorar el rendimiento del programa. Este tipo de mejoras cuando es posible usarlas disminuyen sustancialmente el tiempo de ejecución, ya que en gran medida el consumo total de CPU está en la manipulación de matrices, pero esto no hace paralelo al método de subestructuración secuencial. 111 8.4. Método de Subestructuración en Paralelo La computación en paralelo es una técnica que nos permite distribuir una gran carga computacional entre muchos procesadores. Y es bien sabido que una de las mayores dificultades del procesamiento en paralelo es la coordinación de las actividades de los diferentes procesadores y el intercambio de información entre los mismos. Para hacer una adecuada coordinación de actividades entre los diferentes procesadores, el programa que soporta el método de subestructuración paralelo, usa la misma jerarquía de clases que el método de subestructuración secuencial. Este se desarrolló para usar el esquema maestro-esclavo, de forma tal que el nodo maestro mediante la agregación de un objeto de la clase de Geometría genere la descomposición gruesa del dominio y los nodos esclavos creen un conjunto de objetos FEM2D Rectángulos para que en estos objetos se genere la participación fina y mediante el paso de mensajes (vía MPI) puedan comunicarse los nodos esclavos con el nodo maestro, realizando las siguientes tareas: A) El nodo maestro genera la descomposición gruesa del dominio (supongamos particionado en n×m subdominios) mediante la agregación de un objeto de la clase Geometría, esta geometría es pasada a los nodos esclavos. B) Con esa geometría se construyen los objetos FEM2D Rectángulos (uno por cada subdominio), donde cada subdominio es particionado (supongamos en p × q subdominios). Cada objeto de FEM2D Rectángulos genera la geometría solicitada, regresando las coordenadas de los nodos de frontera del subdominio correspondiente al nodo maestro. C) Con estas coordenadas, el nodo maestro conoce a los nodos de la frontera interior (son estos los que resuelve el método de descomposición de dominio). Las coordenadas de los nodos de la frontera interior se dan a conocer a los objetos FEM2D Rectángulos en los nodos esclavos, transmitiendo sólo aquellos que están en su subdominio. D) Después de conocer los nodos de la frontera interior, cada objeto FEM2D Rectángulos calcula las matrices AΣΣ i , AΣI i , AIΣ i y AII i necesarias para construir el complemento de Schur local S i = AΣΣ i − AΣI i ³ AII i ´ −1 AIΣ i sin realizar comunicación alguna. Al terminar de calcular las matrices se avisa al nodo maestro de la finalización de los cálculos. E) Mediante la comunicación de vectores del tamaño del número de nodos de la frontera interior entre el nodo maestro y los objetos FEM2D Rectángulos, se prepara todo lo necesario para empezar el hmPétodo de gradiente conjugado y resolver el sistema lineal virtual E i=1 S i i uΣ = hPE i=1 bi i . 112 F) Para usar el método de gradiente conjugado, se transmite un vector del tamaño del número de nodos de la frontera interior para que en cada objeto se realicen las operaciones pertinentes y resolver así el sistema algebraico asociado, esta comunicación se realiza de ida y vuelta entre el nodo maestro y los objetos FEM2D Rectángulos tantas veces como iteraciones haga el método. Resolviendo con esto los nodos de la frontera interior uΣi. G) Al término de las iteraciones se pasa la solución uΣi de los nodos de la frontera interior que pertenecen a cada subdominio dentro de cada objeto FEM2D Rectángulos para que se resuelvan los nodos interiores uI i = ³ AII i ´ −1 ³ bI i − AIΣ i uΣi ´ , sin realizar comunicación alguna en el proceso, al concluir se avisa al nodo maestro de ello. I) El nodo maestro mediante un último mensaje avisa que se concluya el programa, terminado así el esquema maestro-esclavo. Del algoritmo descrito anteriormente hay que destacar la sincronía entre el nodo maestro y los objetos FEM2D Rectángulos contenidos en los nodos esclavos, esto es patente en las actividades realizadas en los incisos A, B y C, estas consumen una parte no significativa del tiempo de cálculo. Una parte importante del tiempo de cálculo es consumida en la generación de las matrices locales descritas en el inciso D que se realizan de forma independiente en cada nodo esclavo, esta es muy sensible a la discretización particular del dominio usado en el problema. Los incisos E y F del algoritmo consumen la mayor parte del tiempo total del ejecución al resolver el sistema lineal que dará la solución a los nodos de la frontera interior. La resolución de los nodos interiores planteada en el inciso G consume muy poco tiempo de ejecución, ya que sólo se realiza una serie de cálculos locales previa transmisión del vector que contiene la solución a los nodos de la frontera interior. Este algoritmo es altamente paralelizable ya que los nodos esclavos están la mayor parte del tiempo ocupados y la fracción serial del algoritmo esta principalmente en las actividades que realiza el nodo maestro, estas nunca podrán ser eliminadas del todo pero consumirán menos tiempo del algoritmo conforme se haga más fina la malla en la descomposición del dominio. Para resolver la Ec. (318), usando 513 × 513 nodos (igual al ejemplo de FEM2D Rectángulos secuencial), en la cual se toma una partición rectangular gruesa de 4 × 4 subdominios y cada subdominio se descompone en 128 × 128 subdominios. Usando para los cálculos en un procesador el equipo secuencial y para la parte paralela el cluster heterogéneo resolviendo por el método de gradiente conjugado sin precondicionador, la solución se encontró en 159 iteraciones obteniendo los siguientes valores: 113 Procesadores Tiempo Factor de Aceleración Eficiencia Fracción Serial 1 2943 seg. 2 2505 seg. 1.17485 0.58742 0.70234 3 1295 seg. 2.27258 0.75752 0.16004 4 1007 seg. 2.92254 0.73063 0.12289 5 671 seg. 4.38599 0.87719 0.03499 6 671 seg. 4.38599 0.73099 0.07359 7 497seg. 5.92152 0.84593 0.03035 8 497 seg. 5.92152 0.74019 0.05014 9 359 seg. 8.19777 0.91086 0.01223 10 359 seg. 8.19777 0.81977 0.02442 11 359 seg. 8.19777 0.74525 0.03441 12 359 seg. 8.19777 0.68314 0.04216 13 359 seg. 8.19777 0.63059 0.04881 14 359 seg. 8.19777 0.58555 0.05444 15 359 seg. 8.19777 0.54651 0.05926 16 359 seg. 8.19777 0.51236 0.06344 17 188 seg 15.65425 0.92083 0.00537 Estos resultados pueden ser apreciados mejor de manera gráfica como se muestra a continuación: Figura 15: Métricas de desempeño de 2 a 17 procesadores 114 Primeramente notemos que existe mal balanceo de cargas. La descomposición adecuada del dominio para tener un buen balanceo de cargas se logra cuando se descompone en n×m nodos en la partición gruesa, generándose n ∗m subdominios y si se trabaja con P procesadores (1 para el nodo maestro y P −1 para los nodos esclavos), entonces el balance de cargas adecuado será cuando (P − 1) | (n ∗ m). En nuestro caso se logra un buen balanceo de cargas cuando se tienen 2, 3, 5, 9, 17 procesadores, cuyas métricas de desempeño se muestran a continuación: Figura 16: Métricas de desempeño mostrando sólo cuando las cargas están bien balanceadas (2, 3, 5, 9 y 17 procesadores). En cuanto a las métricas de desempeño, obtenemos que el factor de aceleración en el caso ideal debería de aumentar de forma lineal al aumento del número de procesadores, que en nuestro caso no es lineal pero cumple bien este hecho si están balanceadas las cargas de trabajo. El valor de la eficiencia deberá ser cercano a uno cuando el hardware es usado de manera eficiente, como es en nuestro caso cuando se tiene un procesador por cada subdominio. Y en la fracción serial su valor debiera de tender a cero en el caso ideal, siendo este nuestro caso si están balanceadas las cargas de trabajo, de aquí se puede concluir que la granularidad del problema es gruesa, es decir, no existe una sobrecarga en los procesos de comunicación siendo el cluster una buena 115 herramienta de trabajo para este tipo de problemas. Finalmente las posibles mejoras de eficiencia para el método de subestructuración en paralelo para disminuir el tiempo de ejecución pueden ser: • Balanceo de cargas de trabajo homogéneo. • Al compilar los códigos usar directivas de optimización. • Usar bibliotecas que optimizan las operaciones en el manejo de los elementos de la matriz usando punteros en las matrices densas o bandadas. • El cálculo de las matrices que participan en el complemento de Schur pueden ser obtenidas en paralelo. 8.5. Método de Subestructuración en Paralelo Precondicionado En este método a diferencia del método de subestructuración paralelo, se agrega un precondicionador al método de gradiente conjugado precondicionado visto en la sección (5.3.1) a la hora de resolver el sistema algebraico asociado al método de descomposición de dominio. En este caso por el mal condicionamiento de la matriz, los precondicionadores a posteriori no ofrecen una ventaja real a la hora de solucionar el sistema lineal algebraico. Es por ello que usaremos los precondicionadores a priori vistos en la sección (6.2.1). Estos son más particulares y su construcción depende del proceso que origina el sistema lineal algebraico. Existe una amplia gama de este tipo de precondicionadores, pero son específicos al método de descomposición de dominio usado. Para el método de subestructuración usaremos el derivado de la matriz de rigidez, este no es el precondicionador óptimo para este problema, pero para fines demostrativos nos vasta. La implementación de los métodos a priori, requieren de más trabajo tanto en la face de construcción como en la parte de su aplicación, la gran ventaja de este tipo de precondicionadores es que pueden ser óptimos, es decir, para ese problema en particular el precondicionador encontrado será el mejor precondicionador existente, llegando a disminuir el número de iteraciones hasta en un orden de magnitud. Por ejemplo, al resolver la Ec. (318) usando 513 × 513 nodos en la cual se toma una partición rectangular gruesa de 4×4 subdominios y cada subdominio se descompone en 128 × 128 subdominios. Usando para el cálculo en un procesador el equipo secuencial y para la parte paralela el cluster heterogéneo resolviendo por el método de gradiente conjugado con precondicionador la solución se encontró en 79 iteraciones (una mejora en promedio cercana al 50% con respecto a no usar precondicionador 159 iteraciones) obteniendo los siguientes valores: 116 Procesadores Tiempo Factor de Aceleración Eficiencia Fracción Serial 1 2740 seg. 2 2340 seg. 1.17094 0.58547 0.70802 3 1204 seg. 2.27574 0.75858 0.15912 4 974 seg. 2.81314 0.70328 0.14063 5 626 seg. 4.37699 0.87539 0.03558 6 626 seg. 4.37699 0.72949 0.07416 7 475 seg. 5.76842 0.82406 0.03558 8 475 seg. 5.76842 0.72105 0.05526 9 334 seg. 8.20359 0.91151 0.01213 10 334 seg. 8.20359 0.82035 0.02433 11 334 seg. 8.20359 0.74578 0.03408 12 334 seg. 8.20359 0.68363 0.04207 13 334 seg. 8.20359 0.63104 0.04872 14 334 seg. 8.20359 0.58597 0.05435 15 334 seg. 8.20359 0.54690 0.05917 16 334 seg. 8.20359 0.51272 0.06335 17 173 seg. 15.83815 0.93165 0.00458 Estos resultados pueden ser apreciados mejor de manera gráfica como se muestra a continuación: Figura 17: Métricas de desempeño de 2 a 17 procesadores 117 De las métricas de desempeño, se observa que el factor de aceleración, en el caso ideal debería de aumentar de forma lineal al aumentar el número de procesadores, que en nuestro caso no es lineal pero cumple bien este hecho si está balanceada las cargas de trabajo. En la eficiencia su valor deberá ser cercano a uno cuando el hardware es usado de manera eficiente, como es en nuestro caso cuando se tiene un procesador por cada subdominio. Y en la fracción serial su valor debiera de tender a cero en el caso ideal, siendo este nuestro caso si están balanceadas las cargas de trabajo. En este ejemplo, como en el caso sin precondicionador el mal balanceo de cargas está presente y es cualitativamente igual, para este ejemplo se logra un buen balanceo de cargas cuando se tienen 2, 3, 5, 9, 17 procesadores, cuyas métricas de desempeño se muestran a continuación: Figura 18: Métricas de desempeño mostrando sólo cuando las cargas están bien balanceadas (2, 3, 5, 9 y 17 procesadores). Las mismas mejoras de eficiencia que para el método de subestructuración en paralelo son aplicables a este método, adicionalmente: • Usar el mejor precondicionador a priori disponible para el problema en particular, ya que esto disminuye sustancialmente el número de iteraciones (hasta en un orden de magnitud cuando el precondicionador es óptimo). 118 9. Análisis de Rendimiento y Conclusiones Uno de los grandes retos del área de cómputo científico es poder analizar a priori una serie consideraciones dictadas por factores externos al problema de interés que repercuten directamente en la forma de solucionar el problema, estas consideraciones influirán de manera decisiva en la implementación computacional de la solución numérica. Algunas de estas consideraciones son: • Número de Procesadores Disponibles • Tamaño y Tipo de Partición del Dominio • Tiempo de Ejecución Predeterminado Siendo común que ellas interactúan entre si, de forma tal que normalmente el encargado de la implementación computacional de la solución numérica tiene además de las complicaciones técnicas propias de la solución, el conciliarlas con dichas consideraciones. Esto deja al implementador de la solución numérica con pocos grados de libertad para hacer de la aplicación computacional una herramienta eficiente y flexible que cumpla con los lineamientos establecidos a priori y permita también que esta sea adaptable a futuros cambios de especificaciones -algo común en ciencia e ingeniería-. En este capítulo haremos el análisis de los factores que merman el rendimiento de la aplicación y veremos algunas formas de evitarlo, haremos una detallada descripción de las comunicaciones entre el nodo principal y los nodos esclavos, la afectación en el rendimiento al aumentar el número de subdominios en la descomposición y detallaremos algunas consideraciones generales para aumentar el rendimiento computacional. También daremos las conclusiones generales a este trabajo y veremos las diversas ramificaciones que se pueden hacer en trabajos futuros. 9.1. Análisis de Comunicaciones Para hacer un análisis de las comunicaciones entre el nodo principal y los nodos esclavos en el método de subestructuración es necesario conocer qué se trasmite y su tamaño, es por ello detallaremos en la medida de lo posible las comunicaciones existentes (hay que hacer mención que entre los nodos esclavos no hay comunicación alguna). Tomando la descripción del algoritmo detallado en el método de subestructuración en paralelo visto en la sección (8.4), suponiendo una partición del dominio Ω en n × m y p × q subdominios, las comunicaciones correspondientes a cada inciso son: A) El nodo maestro transmite 4 coordenadas (en dos dimensiones) correspondientes a la delimitación del subdominio. B) 2 ∗ p ∗ q coordenadas transmite cada objeto FEM2D Rectángulos al nodo maestro. 119 C) A lo más n ∗ m ∗ 2 ∗ p ∗ q coordenadas son las de los nodos de la frontera interior, y sólo aquellas correspondientes a cada subdominio son trasmitidas por el nodo maestro a los objetos FEM2D Rectángulos siendo estas a lo más 2 ∗ p ∗ q coordenadas. D) Sólo se envía un aviso de la conclusión del cálculo de las matrices. E) A lo más 2 ∗ p ∗ q coordenadas son transmitidas a los objetos FEM2D Rectángulos desde el nodo maestro y los nodos esclavos transmiten al nodo maestro esa misma cantidad información. F) A lo más 2 ∗ p ∗ q coordenadas son transmitidas a los objetos FEM2D Rectángulos en los nodos esclavos y estos retornan un número igual al nodo maestro por iteración del método de gradiente conjugado. G) A lo más 2 ∗ p ∗ q valores de la solución de la frontera interior son transmitidas a los objetos FEM2D Rectángulos desde el nodo maestro y cada objeto transmite un único aviso de terminación. I) El nodo maestro manda un aviso a cada objeto FEM2D Rectángulos para concluir con el esquema. La transmisión se realiza mediante paso de arreglos de enteros y números de punto flotante que varían de longitud pero siempre son cantidades pequeñas de estos y se transmiten en forma de bloque, por ello las comunicaciones son eficientes. 9.2. Afectación del Rendimiento al Aumentar el Número de Subdominios en la Descomposición Una parte fundamental a considerar es la afectación del rendimiento al aumentar el número de subdominios en descomposición, ya que el complemento de Schur local S i = AΣΣ i − AΣI i ³ AII i ´ −1 AIΣ i involucra el generar las matrices AII i , AΣΣ i , AΣI i , AIΣ i y calcular de alguna forma ³ AII i ´ −1 . Si el número de nodos interiores en el subdominio es grande entonces obtener la matriz anterior será muy costoso computacionalmente, como se ha mostrado en el transcurso de las últimas secciones del capítulo anterior. Al aumentar el número de subdominios en una descomposición particular, se garantiza que las matrices a generar y calcular sean cada vez más pequeñas y fáciles de manejar. Pero hay un límite al aumento del número de subdominio en cuanto a la eficiencia de ejecución, este cuello de botella es generado por el esquema maestroesclavo y es reflejado por un aumento del tiempo de ejecución al aumentar el número de subdominios en una configuración de hardware particular. Esto se debe a que en el esquema maestro-esclavo, el nodo maestro deberá de atender todas las peticiones hechas por cada uno de los nodos esclavos, esto 120 toma especial relevancia cuando todos o casi todos los nodos esclavos compiten por ser atendidos por el nodo maestro. Por ello se recomienda implementar este esquema en un cluster heterogéneo en donde el nodo maestro sea más poderoso computacionalmente que los nodos esclavos. Si a éste esquema se le agrega una red de alta velocidad y de baja latencia, se le permitirá operar al cluster en las mejores condiciones posibles, pero este esquema se verá degradado al aumentar el número de nodos esclavos inexorablemente. Por ello hay que ser cuidadosos en cuanto al número de nodos esclavos que se usan en la implementación en tiempo de ejecución versus el rendimiento general del sistema al aumentar estos. 9.3. Descomposición Óptima para un Equipo Paralelo Dado. Otra cosa por considerar es que normalmente se tiene a disposición un número fijo de procesadores, con los cuales hay que trabajar, así que es necesario encontrar la descomposición adecuada para esta cantidad de procesadores. No es posible hacer un análisis exhaustivo, pero mediante pruebas podemos determinar cual es la mejor descomposición en base al tiempo de ejecución. Para el análisis, consideremos pruebas con 3, 4, 5 y 6 procesadores y veremos cual es la descomposición más adecuada para esta cantidad de procesadores tomando como referencia el resolver la Ec. (318), usando 513 × 513 nodos. Usando para estos cálculos el cluster homogéneo, al resolver por el método de gradiente conjugado precondicionado para cada descomposición se obtuvieron los siguientes resultados: Pa r t i c ión T iem p o en 3 P ro c e sad o re s T iemp o en 4 P ro c esad o re s Tiem p o en 5 Pro ce sad o re s T iemp o en 6 P ro cesado res 2 × 2 y 256 × 256 2576 seg. 2084 seg. 1338 seg. − 4 × 4 y 128 × 128 1324 seg. 1071 seg. 688 seg. 688 seg. 8 × 8 y 64 × 64 779 seg. 630 seg. 405 seg. 405 seg. 16 × 16 y 32 × 32 485 seg. 391 seg. 251 seg. 251 seg. De estas pruebas se observa que el mal balanceo de cargas es reflejado en los tiempos de ejecución, pese a contar con más procesadores no hay una disminución del tiempo de ejecución. Ahora para las mismas descomposiciones, usando el cluster heterogéneo para cada descomposición se obtuvieron los siguientes resultados: Pa r t i c ión T iem p o en 3 P ro c e sad o re s T iemp o en 4 P ro c esad o re s Tiem p o en 5 Pro ce sad o re s T iemp o en 6 P ro cesado res 2 × 2 y 256 × 256 2342 seg. 1895 seg. 1217 seg. − 4 × 4 y 128 × 128 1204 seg. 974 seg. 626 seg. 626 seg. 8 × 8 y 64 × 64 709 seg. 573 seg. 369 seg. 369 seg. 16 × 16 y 32 × 32 441 seg. 356 seg. 229 seg. 229 seg. 121 Primeramente hay que destacar que los nodos esclavos de ambos clusters son comparables en poder de cómputo, pero aquí lo que hace la diferencia es que el nodo maestro del segundo ejemplo tiene mayor rendimiento. Es por ello que al disminuir la fracción serial del problema y atender mejor las comunicaciones que se generan en el esquema maestro-esclavo con todos los objetos FEM2D Rectángulos creados en cada uno de los nodos esclavos mejora sustancialmente el tiempo de ejecución. En ambas pruebas el mal balanceo de cargas es un factor determinante del rendimiento, sin embargo el usó de un cluster en el que el nodo maestro sea más poderoso computacionalmente, hará que se tenga una mejora sustancial en el rendimiento. Para evitar el mal balanceo de cargas se debe de asignar a cada nodo esclavo una cantidad de subdominios igual. La asignación mínima del número de nodos por subdominio queda sujeta a la velocidad de los procesadores involucrados para disminuir en lo posible los tiempos muertos, obteniendo así el máximo rendimiento. La asignación máxima del número de nodos por subdominio a cada nodo esclavo, estará en función de la memoria que consuman las matrices que contienen cada uno de los objetos de FEM2D Rectángulos. La administración de ellos y las comunicaciones no son un factor limitante y por ello se pueden despreciar. Descomposición Fina del Dominio Supongamos ahora que deseamos resolver el problema de una descomposición fina del dominio Ω en 65537 × 65537 nodos, este tipo de problemas surgen cotidianamente en la resolución de sistemas reales y las opciones para implantarlo en un equipo paralelo son viables, existen y son actualmente usadas. Aquí las opciones de partición del dominio son muchas y variadas, y la variante seleccionada dependerán fuertemente de las características del equipo de cómputo paralelo del que se disponga, es decir, si suponemos que una descomposición de 1000 × 1000 nodos en un subdominio consume 1 GB de RAM y el consumo de memoria crece linealmente con el número de nodos, entonces algunas posibles descomposiciones son: Procesadores Descomposición Nodos Subdominio RAM Mínimo 5 2 × 2 y 32768 × 32768 32768 × 32768 ≈33.0 GB 257 16 × 16 y 4096 × 4096 4096 × 4096 ≈4.0 GB 1025 32 × 32 y 2048 × 2048 2048 × 2048 ≈2.0 GB 4097 64 × 64 y 1024 × 1024 1024 × 1024 ≈1.0 GB Notemos que para las primeras particiones, el consumo de RAM es excesivo y en las últimas particiones la cantidad de procesadores en paralelo necesarios es grande (pero ya de uso común en nuestros días). Como en general, contar con equipos paralelos de ese tamaño es en extremo difícil, ¿es posible resolver este tipo de problemas con una cantidad de procesadores menor al número sugerido y donde cada uno de ellos tiene una memoria muy por debajo de lo sugerido?, la respuesta es si. 122 Primero, notemos que al considerar una descomposición del tipo 64 × 64 y 1024 × 1024 subdominios requerimos aproximadamente 1.0 GB de RAM mínimo por nodo, si suponemos que sólo tenemos unos cuantos procesadores con memoria limitada (digamos 2 GB), entonces no es posible tener en memoria de manera conjunta a las matrices generadas por el método. Una de las grandes ventajas de los métodos de descomposición de domino es que los subdominios son en principio independientes entre si y que sólo están acoplados a través de la solución en la interfaz de los subdominios que es desconocida. Como sólo requerimos tener en memoria la información de la frontera interior, es posible bajar a disco duro todas las matrices y datos complementarios (que consumen el 99% de la memoria del objeto FEM2D Rectángulos) generados por cada subdominio que no se requieran en ese instante para la operación del esquema maestroesclavo. Recuperando del disco duro solamente los datos del subdominio a usarse en ese momento (ya que el proceso realizado por el nodo maestro es secuencial) y manteniéndolos en memoria por el tiempo mínimo necesario. Así, es posible resolver un problema de una descomposición fina, usando una cantidad de procesadores fija y con una cantidad de memoria muy limitada por procesador. En un caso extremo, la implementación para resolver un dominio Ω descompuesto en un número de nodos grande es posible implementarla usando sólo dos procesos en un procesador, uno para el proceso maestro y otro para el proceso esclavo, en donde el proceso esclavo construiría las matrices necesarias por cada subdominio y las guardaría en disco duro, recuperándolas conforme el proceso del nodo maestro lo requiera. Nótese que la descomposición del domino Ω estará sujeta a que cada subdominio Ωi sea soportado en memoria conjuntamente con los procesos maestro y esclavo. De esta forma es posible resolver un problema de gran envergadura usando recursos computacionales muy limitados, sacrificando velocidad de procesamiento en aras de poder resolver el problema. Está es una de las grandes ventajas de los métodos de descomposición de dominio con respecto a los otros métodos de discretización tipo diferencias finitas y elemento finito. El ejemplo anterior nos da una buena idea de las limitantes que existen en la resolución de problemas con dominios que tienen una descomposición fina y nos pone de manifiesto las características mínimas necesarias del equipo paralelo para soportar dicha implantación. 9.4. Consideraciones para Aumentar el Rendimiento Algunas consideraciones generales para aumentar el rendimiento son: a) Balanceo de cargas de trabajo homogéneo, si se descompone en n × m subdominios en la partición gruesa se y si se trabaja con P procesadores, entonces el balance de cargas adecuado será cuando (P − 1) | (n ∗ m).), ya que de no hacerlo el rendimiento se ve degradado notoriamente. 123 b) Usar el mejor precondicionador a priori disponible para el problema en particular, ya que esto disminuye sustancialmente el número de iteraciones (hasta en un orden de magnitud cuando el precondicionador es óptimo). c) Usar la biblioteca Lapack++ de licencia GNU que optimiza las operaciones en el manejo de los elementos de la matriz usando punteros en donde se usan matrices densas y bandadas. d) Si se cuenta con un equipo en que cada nodo del cluster tenga más de un procesador, usar bibliotecas (PLAPACK por ejemplo) que permitan paralelizar mediante el uso de procesos con memoria compartida, Pipeline o hilos, las operaciones que involucren a vectores y matrices; como una forma de mejorar el rendimiento del programa. e) Siempre usar al momento de compilar los códigos, directivas de optimización (estas ofrecen mejoras de rendimiento en la ejecución de 30% aproximadamente en las pruebas realizadas), pero existen algunos compiladores con optimizaciones especificas para ciertos procesadores (Intel compiler para 32 y 64 bits) que pueden mejorar aun más este rendimiento (más de 50%). Todas estas mejoras pueden ser mayores si se usa un nodo maestro del mayor poder computacional posible aunado a una red en el cluster de de 1 Gb o mayor y de ser posible de baja latencia, si bien las comunicaciones son pocas, estas pueden generar un cuello de botella sustancial. Por otro lado, hay que hacer notar que es posible hacer usó de múltiples etapas de paralelización que pueden realizarse al momento de implantar el método de descomposición de dominio, estas se pueden implementar conforme se necesite eficiencia adicional, pero implica una cuidadosa planeación al momento de hacer el análisis y diseño de la aplicación y una inversión cuantiosa en tiempo para implantarse en su totalidad, estas etapas se describen a continuación: • El propio método de descomposición de dominio ofrece un tipo particular y eficiente de paralelización al permitir dividir el dominio en múltiples subdominios independientes entre si, interconectados sólo por la frontera interior. • A nivel subdominio otra paralelización es posible, específicamente en el llenado de las matrices involucradas en el método de descomposición de dominio, ya que varias de ellas son independientes. • A nivel de los cálculos, entre matrices y vectores involucrados en el método también se pueden paralelizar de manera muy eficiente. • A nivel del compilador es posible generar el ejecutable usando esquemas de paralelización automático y opciones para eficientizar la ejecución. Por lo anterior es posible usar una serie de estrategias que permitan realizar estas etapas de paralelización de manera cooperativa y aumentar la eficiencia 124 en una factor muy significativo, pero esto implica una programación particular para cada una de las etapas y poder distribuir las tareas paralelas de cada etapa en uno o más procesadores distintos a los de las otras etapas. Notemos finalmente que si se toma el programa de elemento finito y se paraleliza usando sólo directivas de compilación, el aumento del rendimiento es notorio pero este se merma rápidamente al aumentar del número de nodos (esto es debido al aumento en las comunicaciones para mantener y acceder a la matriz del sistema algebraico asociado al método). Pero es aun más notorio cuando el método de descomposición de dominio serial usando las mismas directivas de compilación se paraleliza (sin existir merma al aumentar el número de nodos siempre y cuando las matrices generadas estén el la memoria local del procesador). Esto se debe a que en el método de elemento finito la matriz estará distribuida por todos los nodos usando memoria distribuida, esto es muy costoso en tiempo de cómputo ya que su manipulación requiere de múltiples comunicaciones entre los procesadores, en cambio el el método de descomposición de dominio ya están distribuidas las matrices en los nodos y las operaciones sólo involucran transmisión de un vector entre ellos, minimizando las comunicaciones entre procesadores. Pero aún estos rendimientos son pobres con respecto a los obtenidos al usar el método de descomposición de dominio paralelizado conjuntamente con bibliotecas para manejo de matrices densas y dispersas en equipos con nodos que cuenten con más de un procesador, en donde mediante el uso de memoria compartida se pueden usar el resto de los procesadores dentro del nodo para efectuar en paralelo las operaciones en donde estén involucradas las matrices. 9.5. Conclusiones Generales A lo largo del presente trabajo se ha mostrado que al aplicar métodos de descomposición de dominio conjuntamente con métodos de paralelización es posible resolver una gama más amplia de problemas de ciencias e ingeniería que mediante las técnicas tradicionales del tipo diferencias finitas y elemento finito. La resolución del sistema algebraico asociado es más eficiente cuando se hace uso de precondicionadores a priori conjuntamente con el método de gradiente conjugado precondicionado al implantar la solución por el método de descomposición de dominio. Y haciendo uso del análisis de rendimiento, es posible encontrar la manera de balancear las cargas de trabajo que son generadas por las múltiples discretizaciones que pueden obtenerse para la resolución de un problema particular, minimizando en la medida de lo posible el tiempo de ejecución y adaptándolo a la arquitectura paralela disponible, esto es especialmente útil cuando el sistema a trabajar es de tamaño considerable. Adicionalmente se vieron los alcances y limitaciones de esta metodología, permitiendo tener cotas tanto para conocer las diversas descomposiciones que es posible generar para un número de procesadores fijo, como para conocer el número de procesadores necesarios en la resolución de un problema particular. 125 También se vio una forma de usar los métodos de descomposición de dominio en casos extremos en donde una partición muy fina, genera un problema de gran envergadura y como resolver este usando recursos computacionales muy limitados, sacrificando velocidad de procesamiento en aras de poder resolver el problema. Así, podemos afirmar de manera categórica que conjuntando los métodos de descomposición de dominio, la programación orientada a objetos y esquemas de paralelización que usan el paso de mensajes, es posible construir aplicaciones que coadyuven a la solución de problemas en dos o más dimensiones concomitantes en ciencia e ingeniería, los cuales pueden ser de tamaño considerable. Las aplicaciones desarrolladas bajo este paradigma serán eficientes, flexibles y escalables; a la vez que son abiertas a nuevas tecnologías y desarrollos computacionales y al ser implantados en clusters, permiten una codificación ordenada y robusta, dando con ello una alta eficiencia en la adaptación del código a nuevos requerimientos, como en la ejecución del mismo. De forma tal que esta metodología permite tener a disposición de quien lo requiera, una gama de herramientas flexibles y escalables para coadyuvar de forma eficiente y adaptable a la solución de problemas en medios continuos de forma sistemática. 9.6. Trabajo Futuro Pese a que en este trabajo sólo se mostraron dos métodos de descomposición de dominio, es posible adaptar la metodología usada en este trabajo a muchos otros métodos de descomposición de dominio (como Trefftz-Herrera, FETI, Multigrid, entre otros) donde cada uno de ellos acepta interpoladores de distintos grados. Permitiendo de esta manera tener un grupo de herramientas que pueden ser usadas en múltiples problemas escogiendo la que ofrezca mayores ventajas computacionales, pero las ventajas y desventajas deberán de ser sopesadas al implementarse en un problema particular. Además, a cada método de descomposición de dominio se pueden construir diversos precondicionadores a priori, con el objetivo de obtener un balance entre la complejidad de los diversos precondicionadores (aunado al del método de descomposición de dominio) y el aumento del rendimiento, tomando en cuenta que el precondicionador óptimo para un problema particular puede involucrar mucho trabajo de programación. Finalmente, estos métodos no están constreñidos a problemas elípticos, pueden adaptarse a problemas parabólicos e hiperbólicos, tanto lineales como no lineales, permitiendo así atacar una gran gama de problemas en medios continuos, para más detalles ver [15] y [2]. 126 10. Apéndice En este apéndice se darán algunas definiciones que se usan a lo largo del presente trabajo, así como se detallan algunos resultados generales de algebra lineal y análisis funcional (en espacios reales) que se anuncian sin demostración pero se indica en cada caso la bibliografía correspondiente donde se encuentran estas y el desarrollo en detalle de cada resultado. 10.1. Nociones de Algebra Lineal A continuación detallaremos algunos resultados de álgebra lineal, las demostraciones de los siguientes resultados puede ser consultada en [20]. Definición 31 Sea V un espacio vectorial y sea f (•) : V → R, f es llamada funcional lineal si satisface la condición f (αv + βw) = αf (v) + βf (w) ∀v, w ∈ V y α,β ∈ R. (324) Definición 32 Si V es un espacio vectorial, entonces el conjunto V ∗ de todas las funcionales lineales definidas sobre V es un espacio vectorial llamado espacio dual de V. Teorema 33 Si {v1, ..., vn} es una base para el espacio vectorial V , entonces existe una única base {v∗1, ..., v∗n } del espacio vectorial dual V ∗ llamado la base dual de {v1, ..., vn} con la propiedad de que V ∗ i = δij . Por lo tanto V es isomorfo a V ∗. Definición 34 Sea D ⊂ V un subconjunto del espacio vectorial V. El nulo de D es el conjunto N (D) de todas las funcionales en V ∗ tal que se nulifican en todo el subconjunto D, es decir N (D) = {f ∈ V ∗ | f (v) = 0 ∀v ∈ D} . (325) Teorema 35 Sea V un espacio vectorial y V ∗ el espacio dual de V, entonces a) N (D) es un subespacio de V ∗ b) Si M ⊂ V es un subespacio de dimensión m, V tiene dimensión n, entonces N (M) tiene dimensión n − m en V ∗. Corolario 36 Si V = L ⊕ M (suma directa) entonces V ∗ = N (L) ⊕ N (M) . Teorema 37 Sean V y W espacios lineales, si T (•) : V → W es lineal, entonces el adjunto T ∗ de T es un operador lineal T ∗ : W ∗ → V ∗ definido por T ∗(w∗)(u) = w∗(T u). (326) Teorema 38 Si H es un espacio completo con producto interior, entonces H∗ = H. 127 Definición 39 Si V es un espacio vectorial con producto interior y T (•) : V → V es una transformación lineal, entonces existe una transformación asociada a T llamada la transformación auto-adjunta T ∗ definida como hT u, vi = hu, T ∗vi . (327) Definición 40 Sea V un espacio vectorial sobre los reales. Se dice que una función τ (•, •) : V ×V → R es una forma bilineal sobre V, si para toda x, y, z ∈ V y α, β ∈ R se tiene τ (αx + βy, z) = ατ (x, z) + βτ (y, z) (328) τ (x, αy + βz) = ατ (x, y) + βτ (x, z) . Definición 41 Si τ (•, •) es una forma bilineal sobre V, entonces la función qτ (•) : V → R definida por qτ (x) = τ (x, x) ∀x ∈ V (329) se le llama la forma cuadrática asociada a τ . Notemos que para una forma cuadrática qτ (•) se tiene que qτ (αx) = |α|2 qτ (x) ∀x ∈ V y α ∈ R. 10.2. σ-Algebra y Espacios Medibles A continuación detallaremos algunos resultados conjuntos de espacios σ- algebra, conjuntos de medida cero y funciones medibles, las demostraciones de los siguientes resultados puede ser consultada en [22] y [3]. Definición 42 Una σ-algebra sobre un conjunto Ω es una familia ξ de subconjuntos de Ω que satisface ∅ ∈ ξ Si ψn ∈ ξ entonces ∞S n=1 ψn ∈ ξ Si ψ ∈ ξ entonces ψc ∈ ξ. Definición 43 Si Ω es un espacio topológico, la familia de Borel es el conjunto σ-algebra más pequeño que contiene a los abiertos del conjunto Ω. Definición 44 Una medida μ sobre Ω es una función no negativa real valuada cuyo dominio es una σ-algebra ξ sobre Ω que satisface μ (∅) = 0 y Si {ψn} es una sucesión de conjuntos ajenos de ξ entonces μ à ∞[ n=1 ψn ! = ∞X n=1 μ (ψn) . (330) 128 Teorema 45 Existe una función de medida μ sobre el conjunto de Borel de R llamada la medida de Lebesge que satisface μ ([a, b]) = b − a. Definición 46 Una función f : Ω → R es llamada medible si f −1 (U ) es un conjunto medible para todo abierto U de Ω. Definición 47 Sea E ⊂ Ω un conjunto, se dice que el conjunto E tiene medida cero si μ (E) = 0. Teorema 48 Si α es una medida sobre el espacio X y β es una medida sobre el espacio Y, podemos definir una medida μ sobre X × Y con la propiedad de que μ (A × B) = α (A) β (B) para todo conjunto medible A ∈ X y B ∈ Y. Teorema 49 (Fubini) Si f (x, y) es medible en X × Y entonces Z X×Y f (x, y) dμ = Z X Z Y f (x, y) dβdα = Z Y Z X f (x, y) dαdβ (331) en el sentido de que cualquiera de las integrales existe y son iguales. Teorema 50 Una función f es integrable en el sentido de Riemann en Ω si y sólo si el conjunto de puntos donde f (x) es no continua tiene medida cero. Observación 51 Sean f y g dos funciones definidas en Ω, decimos que f y g son iguales salvo en un conjunto de medida cero si f (x) 6= g(x) sólo en un conjunto de medida cero. Definición 52 Una propiedad P se dice que se satisface en casi todos lados, si existe un conjunto E con μ (E) = 0 tal que la propiedad se satisface en todo punto de Ec. 10.3. Espacios Lp Las definiciones y material adicional puede ser consultada en [12], [18] y [3]. Definición 53 Una función medible f (•) (en el sentido de Lebesgue) es llamada integrable sobre un conjunto medible Ω ⊂ Rn si Z Ω |f | dx < ∞. (332) Definición 54 Sea p un número real con p ≥ 1. Una función u (•) definida sobre Ω ⊂ Rn se dice que pertenece al espacio Lp(Ω) si Z Ω |u(x)|p dx (333) es integrable. 129 Al espacio L2(Ω) se le llama cuadrado integrable. Definición 55 La norma L2(Ω) se define como kukL2(Ω) = μZ Ω |u(x)|2 dx ¶1 2 < ∞ (334) y el producto interior en la norma L2(Ω) como hu, viL2(Ω) = Z Ω u(x)v(x)dx. (335) Definición 56 Si p → ∞, entonces definimos al espacio L∞(Ω) como el espacio de todas las funciones medibles sobre Ω ⊂ Rn que sean acotadas en casi todo Ω (excepto posiblemente sobre un conjunto de medida cero), es decir, L∞(Ω) = {u | |u(x)| ≤ k} (336) definida en casi todo Ω, para algún k ∈ R. 10.4. Distribuciones La teoría de distribuciones es la base para definir a los espacios de Sobolev, ya que permiten definir las derivadas parciales de funciones no continuas, pero esta es coincidente con las derivadas parciales clásica si las funciones son continuas, para mayor referencia de estos resultados ver [12], [18] y [3] Definición 57 Sea Ω ⊂ Rn un dominio, al conjunto de todas las funciones continuas definidas en Ω se denotarán por C0(Ω), o simplemente C(Ω). Definición 58 Sea u una función definida sobre un dominio Ω la cual es no cero sólo en los puntos pertenecientes a un subconjunto propio K ⊂ Ω. Sea K la clausura de K. Entonces K es llamado el soporte de u. Decimos que u tiene soporte compacto sobre Ω si su soporte K es compacto. Al conjunto de funciones continuas con soporte compacto se denota por C0(Ω). Definición 59 Sea Zn + el conjunto de todas las n-duplas de enteros no negativos, un miembro de Zn + se denota usualmente por α ó β (por ejemplo α = (α1, α2, ..., αn). Denotaremos por |α| la suma |α| = α1+α2+...+αn y por Dαu la derivada parcial Dαu = ∂|α|u ∂xα1 1 ∂xα2 2 ...∂xαn n (337) así, si |α| = m, entonces Dαu denota la m-ésima derivada parcial de u. Definición 60 Sea Cm(Ω) el conjunto de todas las funciones Dαu tales que sean funciones continuas con |α| = m. Y C∞(Ω) como el espacio de funciones en el cual todas las derivadas existen y sean continuas en Ω. 130 Definición 61 El espacio D(Ω) será el subconjunto de funciones infinitamente diferenciales con soporte compacto, algunas veces se denota también como C∞0 (Ω). Definición 62 Una distribución sobre un dominio Ω ⊂ Rn es toda funcional lineal continua sobre D(Ω). Definición 63 El espacio de distribuciones es el espacio de todas las funcionales lineales continuas definidas en D(Ω), denotado como D∗(Ω), es decir el espacio dual de D(Ω). Definición 64 Un función f (•) es llamada localmente integrable, si para todo subconjunto compacto K ⊂ Ω se tiene Z K |f (x)| dx < ∞. (338) Ejemplo de una distribución es cualquier función f (•) localmente integrable en Ω. La distribución F asociada a f se puede definir de manera natural como F : D(Ω) → R como hF, φi = Z Ω f φdx (339) con φ ∈ D(Ω). Si el soporte de φ es K ⊂ Ω, entonces |hF, φi| = ¯¯¯¯ Z Ω f φdx ¯¯¯¯ = ¯¯¯¯ Z K f φdx ¯¯¯¯ ≤ sup x∈K |φ| Z Ω |f (x)| dx (340) la integral es finita y hF, φi tiene sentido. Bajo estas circunstancias F es llamada una distribución generada por f. Otro ejemplo de distribuciones es el generado por todas las funciones continuas acotadas, ya que estas son localmente integrables y por lo tanto generan una distribución. Definición 65 Si una distribución es generada por funciones localmente integrables es llamada una distribución regular. Si una distribución no es generada por una función localmente integrable, es llamada distribución singular (ejemplo de esta es la delta de Dirac). Es posible definir de manera natural en producto de una función y una distribución. Específicamente, si Ω ⊂ Rn, u pertenece a C∞(Ω), y si f (•) es una distribución sobre Ω, entonces entenderemos uf por la distribución que satisface h(uf ) , φi = hf, uφi (341) para toda φ ∈ D(Ω). Notemos que la anterior ecuación es una generalización de la identidad Z Ω [u (x) f (x)] φ (x) dx = Z Ω f (x) [u (x) φ (x)] dx (342) la cual se satisface si f es localmente integrable. 131 Derivadas de Distribuciones Funciones como la delta de Dirac y la Heaviside no tienen derivada en el sentido ordinario, sin embargo, si estas funciones son tratadas como distribuciones es posible extender el concepto de derivada de tal forma que abarque a dichas funciones, para ello recordemos que: Teorema 66 La versión clásica del teorema de Green es dada por la identidad Z Ω u ∂v ∂xi dx = Z ∂Ω uvnids − Z Ω v ∂u ∂xi dx (343) que se satisface para todas las funciones u, v en C1(Ω), donde ni es la i-ésima componente de la derivada normal del vector n en la frontera ∂Ω de un dominio Ω. Una versión de la Ec. (343) en una dimensión se obtiene usando la formula de integración por partes, quedando como Z b a uv0dx = [uv] |b a − Z b a vu0dx, u, v ∈ C1 [a, b] (344) como un caso particular de la Ec. (343). Este resultado es fácilmente generalizable a un resultado usando derivadas parciales de orden m de funciones u, v ∈ Cm(Ω) pero remplazamos u por Dαu en la Ec. (343) y con |α| = m, entonces se puede mostrar que: Teorema 67 Otra versión del teorema de Green es dado por Z Ω (Dαu) vdx = (−1)|α| Z Ω uDαvdx + Z ∂Ω h(u, v)ds (345) donde h(u, v) es una expresión que contiene la suma de productos de derivadas de u y v de orden menor que m. Ahora remplazando v en la Ec. (345) por φ perteneciente a D(Ω) y como φ = 0 en la frontera ∂Ω tenemos Z Ω (Dαu) φdx = (−1)|α| Z Ω uDαφdx (346) ya que u es m-veces continuamente diferenciable, esta genera una distribución denotada por u, tal que hu, φi = Z Ω uφdx (347) o, como Dαφ también pertenece a D(Ω), entonces hu, Dαφi = Z Ω uDαφdx (348) 132 además, Dαu es continua, así que es posible generar una distribución regular denotada por Dαu satisfaciendo hDαu, φi = Z Ω (Dαu) φdx (349) entonces la Ec. (346) puede reescribirse como hDαu, φi = (−1)|α| hu, Dαφi , ∀φ ∈ D(Ω). (350) Definición 68 La derivada de cualquier distribución f (•) se define como: La α−ésima derivada parcial distribucional o derivada generalizada de una distribución f es definida por una distribución denotada por Dαf, que satisface hDαf, φi = (−1)|α| hf, Dαφi , ∀φ ∈ D(Ω). Nótese que si f pertenece a Cm ¡ ¯Ω ¢ , entonces la derivada parcial distribucional coincide con la derivada parcial α−ésima para |α| ≤ m. Derivadas Débiles Supóngase que una función u (•) es localmente integrable que genere una distribución, también denotada por u, que satisface hu, φi = Z Ω uφdx (351) para toda φ ∈ D(Ω). Además la distribución u posee derivada distribucional de todos los ordenes, en particular la derivada Dαu es definida por hDαu, φi = (−1)|α| hu, Dαφi , ∀φ ∈ D(Ω). (352) por supuesto Dαu puede o no ser una distribución regular. Si es una distribución regular, entonces es generada por una función localmente integrable tal que hDαu, φi = Z Ω Dαu (x) φ (x) dx (353) y se sigue que la función u y Dαu están relacionadas por Z Ω Dαu (x) φ (x) dx = (−1)|m| Z Ω u (x) Dαφ (x) dx (354) para |α| ≤ m. Definición 69 Llamamos a la función (o más precisamente, a la equivalencia de clases de funciones) Dαu obtenida en la Ec. (354), la α−ésima derivada débil de la función u. Notemos que si u pertenece a Cm ¡ ¯Ω ¢ , entonces la derivada Dαu coincide con la derivada clásica para |α| ≤ m. 133 11. Bibliografía Referencias [1] A. Quarteroni, A. Valli; Domain Decomposition Methods for Partial Differential Equations. Clarendon Press Oxford 1999. [2] A. Toselli, O. Widlund; Domain Decomposition Methods - Algorithms and Theory. Springer, 2005. [3] B. D. Reddy; Introductory Functional Analysis - With Applications to Boundary Value Problems and Finite Elements. Springer 1991. [4] B. F. Smith, P. E. Bj∅rstad, W. D. Gropp; Domain Decomposition, Parallel Multilevel Methods for Elliptic Partial Differential Equations. Cambridge University Press, 1996. [5] B. I. Wohlmuth; Discretization Methods and Iterative Solvers Based on Domain Decomposition. Springer, 2003. [6] I. Foster; Designing and Building Parallel Programs. Addison-Wesley Inc., Argonne National Laboratory, and the NSF, 2004. [7] G. Herrera; Análisis de Alternativas al Método de Gradiente Conjugado para Matrices no Simétricas. Tesis de Licenciatura, Facultad de Ciencias, UNAM, 1989. [8] I. Herrera, M. Díaz; Modelación Matemática de Sistemas Terrestres (Notas de Curso en Preparación). Instituto de Geofísica, (UNAM). [9] I. Herrera;Un Análisis del Método de Gradiente Conjugado. Comunicaciones Técnicas del Instituto de Geofísica, UNAM; Serie Investigación, No. 7, 1988. [10] I. Herrera; Método de Subestructuración (Notas de Curso en Preparación). Instituto de Geofísica, (UNAM). [11] J. II. Bramble, J. E. Pasciak and A. II Schatz. The Construction of Preconditioners for Elliptic Problems by Substructuring. I. Math. Comput., 47, 103-134,1986. [12] J. L. Lions & E. Magenes; Non-Homogeneous Bonduary Value Problems and Applications Vol. I, Springer-Verlag Berlin Heidelber New York 1972. [13] K. Hutter & K. Jöhnk; Continuum Methods of Physical Modeling. Springer- Verlag Berlin Heidelber New York 2004. [14] L. F. Pavarino, A. Toselli; Recent Developments in Domain Decomposition Methods. Springer, 2003. 134 [15] M.B. Allen III, I. Herrera & G. F. Pinder; Numerical Modeling in Science And Engineering. John Wiley & Sons, Inc . 1988. [16] M. Diaz; Desarrollo del Método de Colocación Trefftz-Herrera Aplicación a Problemas de Transporte en las Geociencias. Tesis Doctoral, Instituto de Geofísica, UNAM, 2001. [17] M. Diaz, I. Herrera; Desarrollo de Precondicionadores para los Procedimientos de Descomposición de Dominio. Unidad Teórica C, Posgrado de Ciencias de la Tierra, 22 pags, 1997. [18] P.G. Ciarlet, J. L. Lions; Handbook of Numerical Analysis, Vol. II. North- Holland, 1991. [19] R. L. Burden y J. D. Faires; Análisis Numérico. Math Learning, 7 ed. 2004. [20] S. Friedberg, A. Insel, and L. Spence; Linear Algebra, 4th Edition, Prentice Hall, Inc. 2003. [21] W. Gropp, E. Lusk, A. Skjellem, Using MPI, Portable Parallel Programming Whit the Message Passing Interface. Scientific and Engineering Computation Series, 2ed, 1999. [22] W. Rudin; Principles of Mathematical Analysis. McGraw-Hill International Editions, 1976. [23] X. O. Olivella, C. A. de Sacribar; Mecánica de Medios Continuos para Ingenieros. Ediciones UPC, 2000. [24] Y. Saad; Iterative Methods for Sparse Linear Systems. SIAM, 2 ed. 2000. [25] Y. Skiba; Métodos y Esquemas Numéricos, un Análisis Computacional. UNAM, 2005. 135 Declaro terminado este trabajo sufrido, ideado y llevado a cabo entre los años 2004 al 2006, aún y a pesar de impedimentos tales como: la mala suerte, la desventura, el infortunio, la incomprensión, la gripa, las horas de frío, la tristeza, la desesperanza, el cansancio, el presente, el pasado y mi futuro, el que dirán, la vergüenza, mis propias incapacidades y limitaciones, mis aversiones, mis temores, mis dudas y en fin, todo aquello que pudiera ser tomado por mi, o por cualquiera, como obstáculo en este tiempo de mentiras, verdades, de incredulidad e ignorancia o negación de la existencia real y física de la mala fe. Atentamente Antonio Carrillo Ledesma Manuel Jafet kitim Jiménez casto Conclusiones : esta tesis de investigación es fundamental en el tiempo que vivimos hoy dia se an echo frecuentes los sismos de varias magnitudes los terremotos, tsunamis, y las alarmas que hay dan solo unos cuantos minutos que no es posible evacuar al mayor número de personas la tierra esta en constantes cambios con un mejor modelo de detección de sismos se podría salvar vidas y también estudiar posibles causas y por si fuera poco anticiparse a los sucesos, un proyecto de esta magnitud tendría un impacto impresionante no solo en la sociedad y el mercado laboral sino el el mundo entero. Fuente: file:///C:/Users/jafet/Desktop/AplicacionesDelComputoEnParaleloALaModelacionDeSistemasTerrestres%202.pdf