Processing math: 100%

PESTAÑAS

lunes, 6 de febrero de 2012

Computación en paralelo I

Hola a todos


En esta ocasión voy a hablar muy escuetamente de lo que va a ser el futuro de la programación científica: la computación en paralelo.

Actualmente los procesadores tienden a tener un número de núcleos mayor de uno, tal y como era hace unos años: todos estamos familiarizados con los "Pentium core 2 duo" o los modernos "Pentium I3, I5, etc". Sin embargo, por lo que a mi respecta, la formación universitaria en ingeniería sólo cubre nociones básicas de programación en serie. Es obvio que, desde un punto de vista computacional, si se poseen varios núcleos de procesamiento, el realizar una programación capaz de utilizar simultáneamente la potencia de cálculo de esos cores va a ser mucho más ventajoso.

Sin embargo, hasta hace bien poco desconocía por completo un punto de partida que me pudiera orientar para empezar a caminar en este vasto campo. Pues bien, si queréis una piedra angular de referencia, vuestra palabra mágica es MPI.

El MPI se puede considerar como un anexo a los lenguajes de programación típicos en ciencias (Fortran, C) que permite procesar un código arrancando varios "threads", o hilos de computación. Realmente lo ideal es ejecutar tantos threads como núcleos tengas en tu cluster, pero no tiene por que ser así.

Para realizar un código en paralelo, hay que pasar de una mentalidad de cálculo en serie, y comenzar a pensar como se puede distribuir todas las tareas que tiene que hacer el código: habrá tareas, como una búsqueda de elementos en un fichero o matriz, que se puede "descomponer" en partes que pueden llevar a cabo cada núcleo. En cambio habrá otras que no se puedan descomponer de modo directo, como la escritura y lectura simultánea de un fichero.

Todas estas tareas deben ser coordinadas, esto es, normalmente se programa de forma que uno de los núcleos actúa como maestro y el resto como esclavos: el primero realizará tareas no paralelizables y se encargará de enviar y recibir información que los esclavos necesiten y/o proporcionen. Los cores esclavos se encargarán de realizar la computación respectiva, enviando o recibiendo información vital para el código.

MPI es una alternativa básica a la programación tradicional. Sin embargo, otras aplicaciones están desarrollándose rápidamente, con el afán de simplificar los procesos de comunicación entre cores y la distribución del trabajo: un ejemplo clásico es el OpenMP.


En próximos posts hablaré más en detalle de MPI y OpenMP.


Un saludo

miércoles, 1 de febrero de 2012

Introducción sobre elementos finitos II

Hola a todos

Tras un largo parón navideño (de casi 2 meses), quería volver a retomar mi actividad hablando un poco más de los elementos finitos.

Como recordatorio breve, comentar que el método de los elementos finitos es una técnica matemática que permite obtener el cálculo aproximado de la solución de una ecuación determinada sin llegar a resolverla. Este método permite la implementación de la formulación débil de las ecuaciones diferenciales presentes en los modelos físicos: lo que quiero decir es que este método satisface forma integral de nuestro problema de contorno, y no el problema de contorno en si.

El método se basa, esencialmente, en aproximar nuestra solución por una combinación lineal de ciertas funciones que nosotros propondremos y que se denominan "funciones de base". Esto es, si nuestro problema es el siguiente:

L(u)=f

donde u es nuestra variable y L(u) es un operador aplicado a u, aproximamos la solución de u por \hat{u}:

\hat{u} = \Sigma a_i · N_i

donde N_i son las funciones de forma, que a priori son conocidas.

Existen una variedad muy grande de funciones de forma: normalmente, estas misteriosas funciones que aparecen en todos los libros como la panacea para resolver todos nuestros problemas por arte de magia no son más que el fruto del estudio de lo que se denomina teoría de espacios funcionales.
Este campo se dedica a hallar funciones, a priori desconocidas, que cumplan una serie de propiedades que puedan ser útiles para la resolución de problemas matemáticos: cumplir con condiciones de contorno, con condiciones de integrabilidad, topología de dimensión determinada.

En los libros suelen aparecer familias de funciones más o menos complejas: normamente poseen una forma de spline (funciones polinomiales con n puntos conocidos), si bien todas ellas suelen ser de cuadrado integrable: esta condición asegura que la integral de dicha función no va a ser infinito, lo cual es condición necesaria en la resolución de cualquier problema físico real. Ejemplos clásicos pueden ser los polinomios de Lagrange o los polinomios de Hermite. Sin embargo, otras funciones de base pueden ser utilizados en ciertos problemas, como las series de Fourier.

Por lo tanto, conocida la forma de dichas funciones, el objetivo del método, esencialmente, consistirá en buscar aquellos valores que cumplan que:

L(\hat{u}) -   f  \approx 0

Sabemos que, al no ser u solución exacta, no lo cumple. Como ya dijimos, entonces debe cumplir:

L(\hat{u}) -   f  = r

Con r, función residuo no nulo de la ecuación. La formulación integral del problema entra aquí. Nosotros vamos a imponer que el valor del residuo en promedio a lo largo de toda la función va a ser cero. Para ello, es necesario integrar el valor de r a lo largo del dominio.

\int w(u) r (u) d\Omega = 0

Esta condición se denomina en matemáticas "Método de los residuos ponderados".

La función w(u) se denomina función peso. La pregunta que surge es obvia: ¿Por que utilizamos esta fórmula con una función misteriosa w(u)? ¿Que motivo real existe para utilizar dicha expresión?. Bien, en el fondo la razón que subyace es la de que queremos que nuestra solución sea una buena solución. El método de elementos finitos busca una SOLUCIÓN APROXIMADA a la ecuación diferencial, cuyas soluciones existen en un subespacio concreto. Nosotros no sabemos cual es en concreto (sólo podemos conocer el ERROR que cometemos entre nuestra aproximación y el valor real),  pero si podemos establecer un CRITERIO DE APROXIMACIÓN, esto es, establecer que aproximación puede ser válida para nosotros.
 Esto se realiza matemáticamente definiendo la función residuo r(u) y la función peso w(u) de tal modo que cumpla nuestro criterio de aproximación: por ejemplo, si quiero seguir un criterio de mínima distancia geométrica entre mi solución real y la aproximada, definiré mi función residuo como el cuadrado de la diferencia (utilizaría el método de mínimos cuadrados), o si quiero que mi función valga exactamente cero en una serie de puntos, utilizaré unas funciones peso tipo delta de Dirac en esos puntos (método de colocación). Al establecer este criterio de aproximación estamos acotando matemáticamente un subespacio de soluciones válido dentro del subespacio que forman mis funciones de forma. El método se encargará de "ajustar" tus resultados a las condiciones de tu problema de la forma que hayas elegido.

En el caso del método de Galerkin (uno de los métodos más usados), las funciones peso que se escogen son las propias funciones de forma del subespacio que forma mi solución. Esta condición que se establece con el método de residuos ponderados en Galerkin es totalmente lógica: yo impongo que la mejor solución aproximada de la ecuación diferencial a resolver es su PROYECCIÓN ORTOGONAL sobre mi espacio de soluciones que me he definido con las funciones de forma que he escogido (esto matemáticamente se expresa como que su producto escalar debe ser nulo) . De este modo las funciones de forma se escogen para que sean ortogonales al propio error.

Todo este proceso, que en apariencia es relativamente sencillo, a la hora de realizar una implementación numérica (me refiero a utilizar un ordenador para realizar los cálculos) resulta hartamente más complejo. La razón es la propia naturaleza de los ordenadores: no están diseñados para trabajar con variables continuas. Por ello, lo que para nosotros en la teoría es una función continua, al ordenador hay que introducirlo como una serie discreta de puntos que representen dichas funciones. Esto obliga en la práctica a utilizar una notación matricial para resolver el problema.

En el próximo artículo voy a detallar un poco más del proceso de implementación de elementos finitos en un ordenador.