You are currently browsing the monthly archive for agosto 2012.

A continuación se muestra una imagen con todos los fichero generados de forma automática:

y el código que contienen algunos de estos ficheros:

Anuncios

Cada instanciación de una clase define un objeto: Persona es una clase, Jose es un objeto, una instancia de Persona. La clase es una plantilla que permite crear objetos.

Una clase almacena un estado: un conjunto de atributos que guardan información del objeto; y un comportamiento: un conjunto de métodos que nos permiten trabajar con el objeto.

A continuación, nuestro diagrama de clases con una muestra de atributos, métodos y relaciones:

En el libro “Geometria Diferencial i Relativitat” de J. Girbau, encontramos que una variedad diferenciable de dimensión n y de clase C^k es un par (M,\mathcal{A}) formado por un conjunto M y una estructura diferenciable \mathcal{A} de dimensión n y clase C^k (de aquí en adelante supondremos dimensión n y clase C^k) y tal que M con la topología asociada es Hausdorff y 2AN.

Si no hay posibilidad de confusión se denota la variedad simplemente por M y es usual utilizar la expresión variedad diferenciable cuando es de clase C^\infty.

Llamamos estructura diferenciable a cada una de las clases de equivalencia por la relación C^k\mbox{-compatibles}: dos atlas A_1 y A_2 de clase C^k son C^k\mbox{-compatibles} si A_1 \cup A_2 también es un atlas de clase C^k. Como dos atlas C^k\mbox{-compatibles} definen la misma topología en M, podemos hablar de topología asociada a una estructura diferenciable.

Un atlas de clase C^k es una familia de aplicaciones biyectivas \mathcal{A}=\{ \varphi_i: U_i \rightarrow A_i\}_{i \in I} llamadas cartas locales, donde I es un conjunto de índices, U_i \subset M abierto y A_i \subset \mathbb{R}^n abierto cumpliendo:

  1. \cup_{i \in I} U_i = M
  2. \forall i,j \in I tal que U_i \cap U_j \neq \emptyset, entonces \varphi_i(U_i \cap U_j) y \varphi_j(U_i \cap U_j) son abiertos y las aplicaciones \varphi_j \circ \varphi_i^{-1}:\varphi_i(U_i \cap U_j) \longrightarrow \varphi_j(U_i \cap U_j) son C^k

Si \mathcal{A} es un atlas sobre un conjunto M entonces

\mathcal{B}_{\mathcal{A}} = \{ \varphi_{i}^{-1}(W) : i \in I, W \subset \mathbb{R}^n \mbox{ abierto}\}

es una base de una topología en M.

Una variedad de Riemann es una variedad diferenciable en la que tenemos definida una métrica: un campo tensorial diferenciable g dos veces covariante, un tensor (0,2), simétrico y definido positivo. A g se le llama métrica de Riemann.

Si (M,g) es una variedad de Riemann y x(t) es una curva diferenciable sobre M, entonces la longitud de x(t) entre dos puntos x(a) y x(b) se define como:

\int_a^b \sqrt{g(\dot{x(t)},\dot{x(t)})} dt

Una variedad pseudoriemanniana es aquella en la que la métrica no cumple la propiedad de ser definida positiva y basta con que sea no degenerada.

Una variedad de Lorentz es aquella cuya métrica tiene signatura (n-1, 1). Son especialmente importantes por sus aplicaciones en la física: la teoria de la relatividad general modela el espacio-tiempo mediante una variedad de Lorentz de dimensión n=4 y signatura (3,1).

Cuando modelamos sistemas, ciertos componentes pueden estar relacionados con otros y estas interrelaciones deben ser modeladas. La relación de asociación es el tipo de relación básico entre dos clases. Se representa mediante una linea donde podemos incluir roles, cardinalidades, dirección y restricciones. Suele implementarse como variables de instancia de una clase en otra.

Con respecto a la implementacion, podriamos tener algo como:

//asociación: obj1 y obj2 tienen una relación.
public class Obj1 { void func(Obj2 obj2) {} };

Una relación de agregación es un tipo particular de asociación aparecido después y que permite describir cuando una clase está formada por otras clases. En las agregaciones, el ciclo de vida de una parte es independiente del ciclo de vida del todo, es decir, que si el padre se elimina, no se eliminan sus hijos. En el diagrama aparecen con una punta de flecha en forma de diamante negro en la clase padre. El código podría ser:

//agregación: obj1 posee a obj2 que se lo han prestado: cuando obj1 muere obj2 puede vivir.
public class Obj1 { private Obj2 obj2; Obj1(Obj2 obj2) { this.obj2 = obj2; } }

Una relación de composición permite describir el mismo tipo de relación que las agregaciones pero, por el contrario, la dependencia de sus ciclos de vida es mas fuerte: si el padre se elimina también se eliminan todos sus hijos, aunque una parte puede ser eliminada sin eliminar el todo. En el diagrama aparecen con una punta de flecha en forma de diamante blanco en la clase padre. En código:

//composición: obj1 posee a obj2 y es responsable de su tiempo de vida: cuando obj1 muere obj2 también muere.
public class Obj1 { private Obj2 obj2 = new Obj2(); }

La lista de clases obtenida junto con sus relaciones de generalización:

Imagen

(Diagrama realizado con Visual Paradigm for UML)

El diagrama de clases muestra los bloques con los que construiremos un sistema orientado a objetos. Describen la vida estática del modelo, mostrando que atributos y comportamientos ocurren sin entrar a detallar de que manera se realizan. También muestran las relaciones existentes entre las clases: generalizaciones que reflejan herencias, agregraciones para reflejar composición o uso y asociaciones para mostrar conexiones.

Para representar cada clase utilizamos un rectángulo con tres compartimentos. En el primero va el nombre de la clases, en el segundo sus atributos que pueden ir acompañados de su valor inicial y en el último los métodos junto a sus parámetros. La visibilidad se indica mediante los simbolos – para indicar que un atributo o método es privado, y + para indicar que un atributo o un método es público.

Una relación de generalización indica herencia. Se representa mediante una punta de flecha triangular en la clase padre que generaliza una clase hijo. Un objeto instanciado de la clase hijo heredará de manera implícita los atributos y métodos de su clase padre. En nombre de las clases abstractas va en cursiva.

Volviendo a nuestro diagrama particular, las clases abstractas las hemos indicado temporalmente mediante un borde discontinuo y los colores indican posibles paquetes: verde para clases específicas de SPH, amarillo para numericalMethods, azul para astrophysics y rojo para dataStructures.

Vamos a hacer una lista con todas las clases que creemos que vamos a necesitar para la implementación del método SPH y aplicarlo a simulaciones de problemas astrofísicos. A continuación pensaremos en sus atributos y métodos. Finalmente, interrelacionaremos las clases entre si. Para todo ello, nos apoyaremos en el Visual Paradigm. Empecemos:

  1. Dispondremos de un stage donde ocurre la acción de nuestra simulación. Este contendrá una clase abstracta spatialDecomposition que podremos implementar posteriormente en las clases derivadas bsptree, octree, kdtree, linkedList y que utilizaremos para selección de partículas vecinas incluso con h variable. También necesitamos una clase mesh/grid para ir resolviendo las ecuaciones de campo de Einstein, que constituyen un problema PDE elíptico, e ir obteniendo el campo gravitacional en el que se mueve nuestro fluido. Por tanto, también necesitaremos un ellipticSolver que derivaremos en spectralSolver o finiteDifferenceSolver.
  2. Necesitamos una clases fluid que almacenará toda la información referente a un fluido .
  3. Los fluidos nos permitiran modelar compactObjects como blackHole, neutronStar, strangeStar… y en donde dispondremos de una operacion merge.
  4. Y también la clase particle que permitirá almacenar la información de cada una de las partículas de nuestros fluidos.
  5. La clase kernel nos permitirá evaluar la función kernel. Como ya comentamos en un post, para ofrecer la posibilidad de utilizar diferentes funciones kernel, la clase kernel será abstracta ofreciendo métodos virtuales que implementaremos en las clases derivadas gaussianKernel, quadraticKernel, cubicKernel, quinticKernel, tensorialKernel.
  6. Necesitamos también una clase odeSolver que utilizaremos para resolver las distintas ODE que aparecen al “discretizar” las ecuaciones de la hidrodinámica. Como tenemos diferentes posibilidades, la implementaremos en una clase abstracta cuyos métodos virtuales los implementaran rungeKutta2OdeSolver, rungeKutta4OdeSolver, velocityStormerVerletOdeSolver, leapfrogOdeSolver, modifiedEulerOdeSolver y que nos permitirá utilizar cualquiera de ellos indistintamente.
  7. La clase abstracta eos nos permitirá trabajar con diferentes ecuaciones de estado.
  8. Desde el punto de vista matemático, necesitamos definir una clase tensor y sus casos particulares scalar, vector y matrix.

El problema de Riemann es un caso especial de problema de valor inicial (IVP) en el que la PDE es:

u_t + a u_x = 0

y la condicion inicial (IC):

u(x,0) = u_0(x) = \begin{cases} u_L\mbox{ si } x < 0 \\ u_R \mbox{ si } x > 0 \end{cases}

donde u_L y u_R son dos valores constantes, de manera que tenemos una discontinuidad en x=0. En este caso, la solución es sencilla:

u(x,t) = u_0(x-at) = \begin{cases} u_L \mbox{ si } x-at < 0\\ u_R \mbox{ si } x-at > 0 \end{cases}.

Podemos extender el problema a un conjunto de m PDEs hiperbólicas:

U_t + AU_x = 0 con -\infty < x < \infty, t>0

donde la matriz A es de coeficientes constantes. Al asumir la hiperbolicidad de A, tenemos m valores propios reales \lambda_i y m vectores propios independientes K^{(i)}. En este caso, la IC se escribe como:

U(x,0) = U^{(0)}(x) = \begin{cases} U_L, x<0\\ U_R, x>0 \end{cases}

Los Riemann solvers son métodos numéricos que permiten resolver el problema de Riemann.

La siguiente operación en física:

a+b+c = d+e

solamente la podremos realizar si todas las variables están en las mismas unidades. En física solemos medir longitudes (L), masas (M), tiempos (T), temperaturas (K), etc. La idea es disponer de unas unidades de manera que ciertas constantes universales tomen el valor 1 simplificando así algunas ecuaciones físicas.

Dada la ecuación v=\frac{d}{dt}x podemos escribir la correspondiente ecuación de dimensión [v] = \frac{L}{T} = L T^{-1}. De la misma manera, nos quedaria [a] = L T^{-2} para las aceleraciones.

Sabemos que la velocidad de la luz en el vacío c = 3 \times 10^{10} \frac{cm}{s} en el sistema cegesimal (CGS: cm, g, s) y su ecuación dimensional es [c] = LT^{-1}. Si queremos conseguir c=1, ¿qué unidades de longitud y tiempo necesitamos? Suponiendo que la longitud la continuamos midiendo en cm, aunque lo denotamos ahora con uL, vamos a definir una unidad de tiempo uT de manera que c=1:

3 \times 10^{10} \frac{uL}{s} \times \frac{1}{x} \frac{s}{uT} = 1 \frac{uL}{uT} \Leftrightarrow 1 uT = \frac{1}{3 \times 10^{10}} s

Por lo que en uL y uT hemos conseguido que c=1.

Tenemos que la constante de gravitación universal G = 6.67 \times 10^{-8} \frac{cm^3}{g s^2} en el sistema CGS, con [G] = L^3 M^{-1} T^{-2}. Si hacemos:

6.67 \times 10^{-8} \frac{uL^3}{g s^2} \times \frac{1}{(3 \times 10^{10})^2} \frac{s^2}{uT^2} \times \frac{1}{x} \frac{g}{uM}=1 \frac{uL^3}{uM uT^2}

tenemos que x = \frac{2.2}{3} \times 10^{-28} y nos queda que 1 g = \frac{2.2}{3} \times 10^{-28} uM. Así pués, G=1 en estas nuevas unidades.

Si prodecemos de esta manera hasta conseguir que  c = G = k_B = 1, donde k_B es la constante de Boltzmann, obtenemos las unidades geometrizadas. Si lo hacemos para c =k_B=h=1, donde h es la constante de Planck, obtenemos las unidades naturales.

En C++ podemos definir clases abstractas. Son clases pensadas para definir un comportamiento, es decir, un conjunto de métodos sin implementación. Sus clases derivadas están obligadas a implementar estos métodos, siempre i cuando queramos que sean clases concretas instanciables. Los métodos sin implementación de las clases abstractas reciben el nombre de métodos virtuales.

En el siguiente diagrama de clases UML podemos ver dos pequeños ejemplos de como podríamos utilizar las clases abstractas. En el primer caso, tenemos la clase abstracta kernel para implemetar la función kernel del método SPH. Uno de sus métodos debe devolver un valor a partir de los parámetros d y h. La clase kernel define la signatura del método virtual calcular. Las clases derivadas quadraticKernel, cubicKernel y quinticKernel lo implementan de manera diferente. La clase que lo utiliza tiene como atributo un kernel abstracto, de manera que le podremos asignar cualquiera de las derivaciones concretas. Lo mismo ocurre con la clase abstracta Eos que tiene un método que devuelve la presión. Tenemos dos clases derivadas que lo implementan, cada cual de una manera. La clase cliente con un atributo Eos podrá instanciar cualquiera de las dos clases derivadas.

El siguiente código, generado automáticamente por Visual Paradigm, muestra como queda implementado el UML2 de kernel i cubicKernel en C++:

kernel.h

namespace example {

 class kernel {

 public:
   virtual double calculate(double d, double h) = 0;

 };

}

cubicKernel.h

namespace example {

  class cubicKernel : example::kernel {

    public:
      double calculate(double d, double h);

  };

}

cubicKernel.cpp


#include "cubicKernel.h"

double example::cubicKernel::calculate(double d, double h) {

  //aquí va la implementación concreta del kernel cúbico
  throw "Not yet implemented";

}

LORENE (Langage Objet pour la RElativité NumériquE), tal y como aparece en la página web de la sección Meudon del Observatorio de Paris, es un conjunto de clases C++ pensada para resolver varios problemas que aparecen en relatividad numérica, y de manera mas amplia en astrofísica computacional. Proporciona una serie de herramientas para resolver ecuaciones en derivadas parciales mediante métodos espectrales multidominio.

Las clases de LORENE implementan estructuras básicas, como vectores y matrices, objetos matemáticos mas abstractos, como tensores, y objetos astrofísicos, como estrellas o agujeros negros.

Es un software libre distribuido bajo licencia GNU General Public License.

Utilizando software de ingenieria inversa, podemos obtener el diagrama de clases de LORENE. De esta manera, podemos ver las relaciones existentes entre estas y comprender mejor el software desarrollado.

El siguiente diagrama, obtenido mediante Visual Paradigm for UML, muestra dos subdiagramas de clases, uno relacionado con agujeros negros y otro con estrellas (algunas clases se muestran expandidas con atributos y metodos mientras que otras se muestran sin expandir):

En [Monaghan 2005] se revisan la teoría y las aplicaciones del SPH desde su aparición hasta la actualidad.

El artículo empieza comentando cinco de los puntos fuertes del método SPH:

  1. La advección se trata de manera exacta: si a una partícula se le asigna un color, y se le especifica una velocidad, el transporte del color por la partícula del sistema es exacto.
  2. Cuando trabajamos con mas de un material, cada uno se describe mediante su propio conjunto de partículas de manera que los problemas de interfaz resultan triviales.
  3. Cierra la brecha entre el contínuo y la fragmentación de una forma natural.
  4. La resolución puede depender tanto de la posición como del tiempo.
  5. SPH tiene la ventaja computacional de que la computación ocurre solo donde está realmente la materia, con la consecuente reducción de calculo y almacenamiento.

En [Monaghan 1992] comenta algunos detalles de implementación del método SPH.

La inicialización de un cálculo SPH comienza especificando la masa, posición, velocidad y energía térmica de cada partícula. Si calculamos la densidad mediante la forma alternativa:

\frac{d}{dt}\rho_a = \sum_b m_b v_{ab} \nabla_a W_{ab}

donde v_{ab} = v_a - v_b y que tiene algunas ventajas (evita la oscilación de los bordes del fluido y es computacionalmente mas eficiente el calculo de la tasa de cambio de todas las variables físicas), necesitamos asignar la densidad inicial de cada partícula. También la mezcla de elementos de cada partícula puede ser necesaria.

Para inicializar las partículas, es conveniente situarlas en una malla regular, que puede ser una malla Cartesiana regular con igual espaciado o una malla cúbica centrada en el cuerpo en tres dimensiones (mejor representación de las integrales por sumas, mejor detección de vecinos y mejor distribución de partículas en caso de compresión planar). Además, si el tamaño de celda asociado con una partícula a es \Delta V_a entonces m_a se toma como \rho \Delta V_a.

Si utilizamos la misma h para cada partícula, entonces basta trabajar con listas enlazadas. Si el kernel esta basado en splines, las celdas enlazadas por la lista deben ser de 2h de ancho de manera que solo las celdas vecinas pueden contribuir en las partículas en una celda dada.

Si cada partícula tiene su propia h entonces el calculo de las sumas del SPH pueden formar parte de un calculo en arbol que es el método natural para calcular las fuerzas autogravitantes de un conjunto de partículas. Este codigo en arbol puede ser vectorizado.

Trabajando en GRMHD nos encontramos el siguiente problema de Cauchy:

  1. Ecuaciones de Einstein: (\gamma,K)_{t=0} \longrightarrow (\gamma,K)_t
  2. Hidrodinámica: (\rho,q,\tau)_{t=0} \longrightarrow (\rho,q,\tau)_t
  3. Electromagetismo: (E,B)_{t=0} \longrightarrow (E,B)_t

donde \gamma es la métrica y K la curvatura extrinseca; \rho es la densidad de masa, q la densidad de momento, \tau la densidad de energia; E es el campo de fuerza eléctrica y B la inducción magnética.

En [Rosswog 2009] comenta algunos métodos para ODEs, ya que una de las características relevantes de SPH es que transforma las PDEs de la hidrodinámica en un sistema de ODEs y hacer una simulación en SPH no es mas que resolver dicho sistema.

Las ODEs en astrofísica suelen ser de segundo orden pero podemos transformarlas en un sistema de ODEs de primer orden introduciendo como nuevas variables a las velocidades. De esta manera obtenemos n ecuaciones acopladas de la forma:

\frac{d}{dt}y_i = f_i(t,y_1,\ldots,y_n) con i = 1,\ldots,n

Tenemos, en este caso, dos tipos de problemas: los problemas de contorno y los problemas de valor inicial. Este último caso, que es el que nos interesa, se trata de determinar el valor de y_i^n en t^n a partir de la condición inicial y_i^0 en t^0.

El método de Euler es el método de integración pedagógico. La solución avanza de t^n a t^{n+1} \equiv t^n + \Delta t según:

\vec{y}^{n+1} = \vec{y}^n + \Delta t \vec{f}(t^n,\vec{y}^n)

donde \vec{y}^{n+1} = \vec{y}(t^{n+1}) y \Delta t es el paso de tiempo escogido. Es un método explícito de primer orden y es asimétrico en tiempo, pues la derivada es correcta al principio del paso de tiempo pero no al final. De ahí el nombre de Euler hacia adelante. En el caso de posiciones, velocidades y aceleraciones, tenemos:

\vec{x}^{n+1} = \vec{x}^n + \Delta t \vec{v}^n

\vec{v}^{n+1} = \vec{v}^n + \Delta t \vec{a}^n

La posibilidad contraria es:

\vec{y}^{n+1} = \vec{y}^n + \Delta t \vec{f}(t^{n+1},\vec{y}^{n+1})

llamado el método de Euler hacia atrás, que es un método implícito. Todos los siguientes métodos que comentaremos son explícitos.

El siguiente método es el Euler modificado que utiliza el valor de la derivada al principio y al final del paso de tiempo:

\vec{y}^{n+1} = \vec{y}^n + \Delta t \frac{\vec{f}^n + \vec{f}^{n+1,p}}{2}

que es un método de segundo orden. En el caso de posiciones, velocidades y aceleraciones tenemos:

\vec{x}^{n+1} = \vec{x}^n + \Delta t \frac{\vec{v}^n + \vec{v}^{n+1,p}}{2}

\vec{v}^{n+1} = \vec{v}^n + \Delta t \frac{\vec{a}^n + \vec{a}^{n+1,p}}{2}

En el método de Runge-Kutta de segundo orden tomamos la derivada en el punto medio de los intervalos de tiempo tomados:

  1. Tomamos como paso \frac{\Delta t}{2}, calculamos \vec{y}^{mid} = \vec{y}^n + \frac{1}{2} \Delta t \vec{f}(t^n,\vec{y}^n) y evaluamos la derivada en ese punto \vec{f}^{mid} = \vec{f}(t^{\frac{n+1}{2}},\vec{y}^{mid})
  2. Tomamos el paso entero con esta derivada: \vec{y}^{n+1} = \vec{y}^n + \Delta t \vec{f}^{mid}

En el caso de posiciones, velocidades y aceleraciones tenemos:

\vec{x}^{n+1} = \vec{x}^n + \Delta t \vec{v}^{mid}

\vec{v}^{n+1} = \vec{v}^n + \Delta t \vec{a}^{mid}

Runge-Kutta de cuarto orden:

\vec{k_1} = \Delta t \vec{f}(t^n,\vec{y}^n)

\vec{k_2} = \Delta t \vec{f}(t^n + \frac{\Delta t}{2},\vec{y}^n + \frac{\vec{k}_1}{2})

\vec{k_3} = \Delta t \vec{f}(t^n + \frac{\Delta t}{2},\vec{y}^n + \frac{\vec{k}_2}{2})

\vec{k_4} = \Delta t \vec{f}(t^n + \Delta t,\vec{y}^n + \vec{k}_3)

\vec{y}^{n+1} = \vec{y}^n + \frac{1}{6}(\vec{k}_1+ 2 \vec{k}_2 + 2 \vec{k}_3 + \vec{k}_4)

Velocity Stormer-Verlet:

En el caso de posiciones, velocidades y aceleraciones tenemos:

\vec{x}^{n+1} = \vec{x}^n + \vec{v}^n \Delta t + \frac{1}{2} \vec{a}^n \Delta t^2

\vec{v}^{n+1} = \vec{v}^n + \Delta t \frac{\vec{a}^n + \vec{a}^{n+1}}{2}

Dado que la conservación exacta es en nuestro puede ser mas importante  que el orden alto, aconseja el último método.

Para un observador inicial, si denotamos con \vec{E} al campo eléctrico, \vec{B} al campo magnético, \rho a la densidad de carga y \vec{J} a la densidad de corriente, entonces tenemos las ecuaciones de Maxwell:

\nabla \cdot \vec{E} = \rho

\nabla \times \vec{E} + \vec{B}_t = 0

\nabla \cdot \vec{B} = 0

\nabla \times \vec{B} - \vec{E}_t = \vec{J}

y la ecuación de continuidad o de conservación de carga:

\rho_t + \nabla \cdot \vec{J} = 0

Si el observador inercial se encuentra en el espacio-tiempo de Minkowski, las ecuaciones de Maxwell se expresan como dos ecuaciones de ligadura:

\nabla \cdot \vec{E} = \rho

\nabla \cdot \vec{B} = 0

y seis ecuaciones de evolución:

\vec{E}_t = \nabla \times \vec{B} - \vec{J}

\vec{B}_t = -\nabla \times \vec{E}

Si en un instante t=t_0 se cumplen las ecuaciones de ligadura y si la carga eléctrica se conserva en un entorno de t=t_0,

\rho_t + \nabla \cdot \vec{J} = 0,

entonces las ecuaciones de ligadura se cumplen en ese entorno (como consecuencia de las ecuaciones de evolución).

En [Rosswog 2009] tenemos las ecuaciones de la hidrodinámica en forma Lagrangiana discretizadas y en su forma mas básica. Las partículas avanzaran en el tiempo siguiendo las siguientes ecuaciones:

Para empezar, no hay necesidad de resolver la ecuación de continuidad ya que la masa de las partículas permanece fija. Podemos obtener las densidades mediante:

\rho_a = \sum_b m_b W_{ab}

La ecuación del momento queda:

\frac{d}{dt}\vec{v}_a = - \sum_b m_b (\frac{P_a}{\rho_a^2} + \frac{P_b}{\rho_b^2} + \Pi_{ab} ) \nabla_a W_{ab}

La ecuación de evolución para la energía interna específica puede escribirse como:

\frac{d}{dt} u_a = \sum_b m_b (\frac{P_a}{\rho_a^2} + \frac{1}{2}\Pi_{ab}) \vec{v}_{ab} \nabla_b W_{ab}

Rosswog las llama “vanilla ice” SPH.

En las Jornadas sobre los problemas del milenio celebradas en Barcelona del 1 al 3 de junio de 2011, Diego Córdoba dió una charla sobre el problema Clay de las ecuaciones de Navier-Stokes para la que escribió estas notas.

Un fluido es ideal si es incompresible, homogéneo  y perfecto.

La idea del problema consiste en determinar si:

…un fluido incompresible con energía finita puede desarrollar singularidades en tiempo finitos.

Mas formalmente, si consideramos un fluido viscoso (\nu > 0), homogéneo (\rho = 1) e incompresible (\nabla \cdot u = 0), tenemos:

u_t + u \cdot \nabla u = - \nabla p + \nu \Delta u + f con \nu >0, x \in \mathbb{R}^3, t \geq 0

\nabla \cdot u = 0

u(x,0) = u_0

donde cada partícula en el tiempo t está en la posición x = (x_1,x_2,x_3) del dominio que ocupa el fluido \Omega \subset \mathbb{R}^3, u(x,t) = (u_1(x,t), u_2(x,t), u_3(x,t)) es el campo de velocidades, p = p(x,t) son las presiones en el seno del fluido y \rho = \rho(x,t) es la densidad. Además, \nu = cte \geq 0 es la viscosidad y f la fuerza externa.

La fuerza exterior  debe verificar:

|\partial_x^\alpha\partial_t^m f| \leq C_{\alpha,k,m}(1+|x|+t)^{-k} para todo \alpha, k, m > 0

y el dato inicial las siguientes condiciones de regularidad:

|\partial_x^\alpha u_{0i}| \leq C_{\alpha,k}(1+|x|)^{-k} para todo \alpha, k > 0

Las soluciones (u,p) admisibles para x \in \mathbb{R}^3 estan o en C^{\infty}(\mathbb{R}^3 \times [0,\infty)) con decaimiento en el infinito de la presión y de energía finita, es decir, \int_{\mathbb{R}^3} |u|^2 dx < \infty para todo t, o estan en C^{\infty}(\mathbb{T}^3 \times [0,\infty)) periódicas y la presión de media cero.

Sea u_0 satisfaciendo las condiciones de regularidad. El problema del Instituto Clay consiste, entonces, en determinar si siempre existen soluciones admisibles para u_0 o si existe algún caso en que no.

Eclipse es un IDE multiplataforma Open Source desarrollado por la Eclipse Foundation. Utiliza módulos, plug-in, para proporcionar diferentes funcionalidades, ya que se basa en la filosofía RCP (Rich Client Platform), de manera que tenemos una plataforma ligera en la que el usuario va añadiendo lo que realmente necesita.

El plug-in CDT (C/C++ Development Toolkit) permite desarrollar en Eclipse proyectos basados en estos lenguajes de programación.

Otro plug-in que encontramos es el MDT (Model Development Tools), el cual permite trabajar, entre otras cosas, con UML2 con el que especificar y diseñar software basado en el paradigma de la orientación a objetos.

En [Monaghan 1992] comenta el caso del método SPH en relatividad especial.

Para empezar asumimos que el fluido está constituido por bariones, por lo que el tensor de energia-momento es:

T^{\mu \nu} = (n m_0 c^2 + n \tau + P) U^\mu U^\nu + P g^{\mu \nu}

donde los indices griegos van de 0 a 3 y los coeficientes de la métrica se definen

g_{00} = -1 y g_{ij} = 1

En estas ecuaciones, n representa la densidad de bariones, P es la presión, \tau la energía térmica, c la velocidad del sonido, U^\nu la 4-velocidad con U_\nu U^\nu = -1 y m_0 la masa.

Las ecuaciones del momento se siguen de:

\frac{\partial}{\partial x^\nu} T^{i \nu} = 0

que es:

\frac{d}{dt} q = - \frac{1}{N} \nabla P

y en forma SPH queda:

\frac{d}{dt}q_a = -\sum_b \nu_b (\frac{P_a}{n_a^2} + \frac{P_b}{N_b^2}) \nabla_a W_{ab}

donde \nu_b es el número de bariones asociados a la partícula b.

Y la de la energía se sigue de:

\frac{\partial}{\partial x^j} T^{0j} = 0

que es:

\frac{d}{dt} \epsilon = - \frac{1}{N} \nabla \cdot (Pv)

y en foma SPH queda:

\frac{d}{dt} \epsilon_a = -\sum_b m_b (\frac{P_a v_a}{N_a^2} + \frac{P_b v_b}{N_b^2}) \nabla W_{ab}

Según Monaghan en [Monaghan, 1992], uno de los padres del método, buscaban un método con el que fuera fácil trabajar y diera buenos resultados y encontraron eso y mucho mas.

Además de no  necesitar malla para calcular las derivadas, puesto que se calculan de manera analítica a partir de una fórmula de interpolación, las ecuaciones de conservación del momento y la energía pasan a ser un conjunto de ecuaciones diferenciales ordinarias fáciles de entender.

En esencia, el método SPH es un método de interpolación que permite expresar una función a partir de los valores en un conjunto desordenado de puntos a los que llamamos particulas.

La idea es que aproximamos cualquier función A(r) de la siguiente manera:

A_I(r) = \int A(r') W(r-r',h) dr'

donde r es el vector posición, W es una función de pesos o kernel (con las propiedades ya comentadas en otros posts) y h es la longitud de suavizado. Este tipo de aproximaciones reciben el nombre de interpolación integral.

Para trabajar de manera numérica, la interpolador integral se aproxima por:

A_S(r) = \sum_b m_b \frac{A_b}{\rho_b} W(r-r_b,h)

donde el índice del sumatorio b identifica a cada partícula y el sumatorio es sobre todas las partículas. La partícula b tiene masa m_b, posición r_b, densidad \rho_b y velocidad v_b. El valor de la función A en r_b es A_b.

En el artículo [Rosswog 2009], Stephan Rosswog hace un repaso detallado del método SPH centrandose especialmente en sus aplicaciones en astrofísica. Repasamos el apartado que hace referencia a las ecuaciones de la hidrodinámica en forma Lagrangiana.

A diferencia de los metodos basados en malla, que son Eulerianos, es decir, métodos donde  describimos el fluido desde un punto fijo del espacio, el Smoothed Particle Hydrodynamics es totalmente Lagrangiano, por lo que describimos el fluido desde un sistema de coordenadas fijado en una particula del fluido en movimiento.

La derivada Lagrangiana o sustancial respecto del tiempo, \frac{d}{dt}, se relaciona con la derivada Euleriana respecto al tiempo, \frac{\partial}{\partial t} de la siguiente manera:

\frac{d}{dt} = \frac{dx^i}{dt} \frac{\partial}{\partial x^i} + \frac{\partial}{\partial t} = \vec{v} \cdot \nabla + \frac{\partial}{\partial t}

Aplicando esta relación a las ecuaciones en forma Euleriana, las ecuaciones de la hidrodinámica en forma Lagrangiana quedan:

  1. Ecuación de continuidad: \frac{d}{dt} \rho = - \rho \nabla \cdot \vec{v}
  2. Ecuacion del momento: \frac{d}{dt} \vec{v} = -\frac{\nabla P}{\rho} + \vec{f}
  3. Ecuación de la energía: \frac{d}{dt}u = \frac{P}{\rho^2} \frac{d}{dt} \rho = - \frac{P}{\rho} \nabla \cdot \vec{v}
  4. Ecuación de estado, que describe la termodinámica del fluido estelar: P = (\gamma -1) \cdot \rho \cdot \epsilon (ecuación del gas ideal)

La hipótesis de Rieman dice que todos los ceros no triviales de la función zeta de Riemann:

\zeta(z) = \sum_{n=1}^{\infty} \frac{1}{n^z}

tienen parte real \frac{1}{2}.

En el entretenido libro “La música de los números primos” de Marcus de Sautoy, éste comenta la posibilidad de que la teoría de números y la física estén mas estrechamente relacionadas de lo que pensamos:

Los físicos creen que la razón por la que los ceros de Riemann deben situarse todos sobre la recta es que terminarán por ser las frecuencias de un tambor matemático. A un cero que se situara fuera de la recta le correspondería una frecuencia imaginaria prohibida por la teoría.

Y para afianzar esta idea, comenta un problema clásico de hidrodinámica resuelto por Bernhard Riemann argumentando de la misma manera:

El problema se refiere a una esfera de fluido en rotación que se mantiene unida gracias a interacciones gravitacionales recíprocas entre las partículas que la componen. Una estrella, por ejemplo, es una enorme bola de gas giratorio que se mantiene unido por su propia gravedad. La cuestión es: ¿qué sucederá con la bola si se le da una patada?¿Se limitará a temblar ligeramente o se desitegrará? Para responder a estas preguntas es necesario determinar si ciertos números imaginarios determinados están o no alineados. Si lo están, la esfera de fluido en rotación quedará intacta.

Los Nachlass son un conjunto de notas inéditas de Riemann que están en la biblioteca de Gotinga. Cuenta el libro que cuando el físico Jon Keating pidió dos partes de los Nachlass, una correspondiente a sus intentos de demostración de la hipotesis de Riemann y otra a sus estudios en hidrodinámica, el bibliotecario le entrego un único grupo de documentos. De ahí que Marcus escriba:

Una vez más, los Nachlass revelaban hasta qué punto Riemann se adelantó a su tiempo. Es imposible que no fuera consciente del significado que implicaba su solución al problema de dinámica de fluidos. Su método había demostrado por qué cietos números imaginarios que aparecían en su análisis de la esfera de fluido se colocaban en línea recta; y al mismo tiempo -y en los mismos folios- estaba intentando demostrar por qué los ceros de su paisaje zeta se situaban todos sobre la misma línea.

 

Las heramientas CASE (Computer Aided Software Engineering) pretenden aplicar los métodos de la ingenieria a la producción de software y surgen en respuesta a la crisis del software de finales de los sesenta. Cubren todas las etapas del ciclo de vida del software y básicamente pretende generar software de calidad minimizando tiempo y costes tanto en el desarrollo como en el mantenimiento del mismo.

El lenguaje UML (Unified Modeling Language) es un estandar que permite especificar y diseñar software basado en el paradigma de la orientación a objetos.

En este enlace podemos ver  diferentes herramientas CASE que soportan UML y que permiten  generar código automáticamente en diferentes lenguajes de programación.

Existen diferentes posibilidades a la hora de definir una función kernel:

  1. Gaussiana [Gingold & Monaghan, 1977]:W(r,h) = \alpha_D \cdot e^{-q^2} con 0 \leq q \leq 2 donde q=\frac{r}{h}, r es la distancia entre dos partículas determinadas y \alpha_D, el factor dimensional, que es \frac{1}{\pi h^2} en dos dimensiones y \frac{1}{\pi^{\frac{3}{2}} h^3} en tres.
  2. Cuadrática [Gingold & Monaghan, 1977]: \alpha_D (\frac{3}{16} q^2 - \frac{3}{4} q + \frac{3}{4}) con 0 \leq q \leq 2 y donde \alpha_D es \frac{2}{\pi h^2} en 2D y \frac{5}{4 \pi h^3} en 3D.
  3. B-spline cúbico [Monaghan & Lattancio, 1985]W(r,h) = \alpha_D \begin{cases} 1 - \frac{3}{2} q^2 + \frac{3}{4} q^3 ,& \mbox{if } 0 \leq q \leq 1 \\ \frac{1}{4}(2-q)^3 ,& \mbox{if } 1 < q \leq 2 \\ 0 ,& \mbox{if } q > 2\end{cases} con \alpha_D = \frac{10}{7 \pi h^2} en dos dimensiones y \frac{1}{\pi h^3} en 3 dimensiones.
  4. Quíntica: W(r,h) = \alpha_D (1-\frac{q}{2}^4)(2q + 1) con 0 \leq q \leq 2 y \alpha_D = \frac{7}{4 \pi h^2} en 2D y \frac{7}{8 \pi h^3} en 3D.

Pensando en posibles implementaciones, parece lógico utilizar una classe abstracta Kernel de la que heredaran las anteriores implementando sus métodos, de manera que, las posibles clases que llamen a la misma puedan hacerlo siempre de la misma manera independientement de cual estemos utilizando.

El siguiente código muestra la definición de una clase en el lenguaje C++. En el podemos ver los elementos típicos de las mismas.

En el fichero Complex.h, tendriamos la representación del tipo:


#ifndef ComplexH
#define ComplexH

class Complex {

  private:

    //atributos
    double re;
    double im;

  public:

    //constructores
    Complex();
    Complex(double re, double im);

    //metodos
    double get_re();
    double get_im();
    void set_re(double re);
    void set_im(double im);
    double get_mod();
    double get_arg();

    //destructor
    ~Complex();

};

Mientras que en el fichero Complex.cpp tendriamos la implementación del mismo:

#include "Complex.h"

//constructores
Complex::Complex() {
  re = 0;
  im = 0;
}

Complex::Complex(double re, double im) {
  this.re = re;
  this.im = im;
}

//destructor
Complex::~Complex() {
}

//metodos
double Complex::get_re() {
  return this.re;
}

double Complex::get_im() {
  return this.im;
}

void Complex::set_re(double re) {
  this.re = re;
}

void Complex::set_im(double im) {
  this.im = im;
}

double Complex::get_mod() {
  return ...;
}

double Complex::get_arg(){
  return ...;
}

};
agosto 2012
L M X J V S D
    Sep »
 12345
6789101112
13141516171819
20212223242526
2728293031