FireStats error : FireStats: Unknown commit strategy

agosto 2012

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:

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 $latex n$ y de clase $latex C^k$ es un par $latex (M,mathcal{A})$ formado por un conjunto $latex M$ y una estructura diferenciable $latex mathcal{A}$ de dimensión $latex n$ y clase $latex C^k$ (de aquí en adelante supondremos dimensión $latex n$ y clase $latex C^k$) y tal que $latex M$ con la topología asociada es Hausdorff y 2AN.

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

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

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

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

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

$latex 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 $latex M$.

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

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

$latex 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 $latex (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 $latex n=4$ y signatura $latex (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 $latex 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:

$latex u_t + a u_x = 0$

y la condicion inicial (IC):

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

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

$latex 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 $latex m$ PDEs hiperbólicas:

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

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

$latex 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:

$latex 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 ($latex L$), masas ($latex M$), tiempos ($latex T$), temperaturas ($latex K$), etc. La idea es disponer de unas unidades de manera que ciertas constantes universales tomen el valor $latex 1$ simplificando así algunas ecuaciones físicas.

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

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

$latex 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 $latex uL$ y $latex uT$ hemos conseguido que $latex c=1$.

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

$latex 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 $latex x = frac{2.2}{3} times 10^{-28}$ y nos queda que $latex 1 g = frac{2.2}{3} times 10^{-28} uM$. Así pués, $latex G=1$ en estas nuevas unidades.

Si prodecemos de esta manera hasta conseguir que  $latex c = G = k_B = 1$, donde $latex k_B$ es la constante de Boltzmann, obtenemos las unidades geometrizadas. Si lo hacemos para $latex c =k_B=h=1$, donde $latex 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 $latex d$ y $latex 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

[sourcecode languaje=”cpp”]
namespace example {

class kernel {

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

};

}
[/sourcecode]

cubicKernel.h

[sourcecode languaje=”cpp”]
namespace example {

class cubicKernel : example::kernel {

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

};

}
[/sourcecode]

cubicKernel.cpp

[sourcecode languaje=”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";

}
[/sourcecode]

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:

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

donde $latex 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 $latex a$ es $latex Delta V_a$ entonces $latex m_a$ se toma como $latex rho Delta V_a$.

Si utilizamos la misma $latex 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 $latex 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 $latex 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: $latex (gamma,K)_{t=0} longrightarrow (gamma,K)_t$
  2. Hidrodinámica: $latex (rho,q,tau)_{t=0} longrightarrow (rho,q,tau)_t$
  3. Electromagetismo: $latex (E,B)_{t=0} longrightarrow (E,B)_t$

donde $latex gamma$ es la métrica y $latex K$ la curvatura extrinseca; $latex rho$ es la densidad de masa, $latex q$ la densidad de momento, $latex tau$ la densidad de energia; $latex E$ es el campo de fuerza eléctrica y $latex 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 $latex n$ ecuaciones acopladas de la forma:

$latex frac{d}{dt}y_i = f_i(t,y_1,ldots,y_n)$ con $latex 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 $latex y_i^n$ en $latex t^n$ a partir de la condición inicial $latex y_i^0$ en $latex t^0$.

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

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

donde $latex vec{y}^{n+1} = vec{y}(t^{n+1})$ y $latex 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:

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

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

La posibilidad contraria es:

$latex 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:

$latex 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:

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

$latex 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 $latex frac{Delta t}{2}$, calculamos $latex 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 $latex vec{f}^{mid} = vec{f}(t^{frac{n+1}{2}},vec{y}^{mid})$
  2. Tomamos el paso entero con esta derivada: $latex vec{y}^{n+1} = vec{y}^n + Delta t vec{f}^{mid}$

En el caso de posiciones, velocidades y aceleraciones tenemos:

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

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

Runge-Kutta de cuarto orden:

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

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

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

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

$latex 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:

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

$latex 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 $latex vec{E}$ al campo eléctrico, $latex vec{B}$ al campo magnético, $latex rho$ a la densidad de carga y $latex vec{J}$ a la densidad de corriente, entonces tenemos las ecuaciones de Maxwell:

$latex nabla cdot vec{E} = rho$

$latex nabla times vec{E} + vec{B}_t = 0$

$latex nabla cdot vec{B} = 0$

$latex nabla times vec{B} – vec{E}_t = vec{J}$

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

$latex 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:

$latex nabla cdot vec{E} = rho$

$latex nabla cdot vec{B} = 0$

y seis ecuaciones de evolución:

$latex vec{E}_t = nabla times vec{B} – vec{J}$

$latex vec{B}_t = -nabla times vec{E}$

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

$latex 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).

Tags: ,

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:

$latex rho_a = sum_b m_b W_{ab}$

La ecuación del momento queda:

$latex 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:

$latex 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.

Tags: , , , , , ,

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 ($latex nu > 0$), homogéneo ($latex rho = 1$) e incompresible ($latex nabla cdot u = 0$), tenemos:

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

$latex nabla cdot u = 0$

$latex u(x,0) = u_0$

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

La fuerza exterior  debe verificar:

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

y el dato inicial las siguientes condiciones de regularidad:

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

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

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

Tags: , ,

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.

Tags: , , , ,

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:

$latex 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 $latex 0$ a $latex 3$ y los coeficientes de la métrica se definen

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

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

Las ecuaciones del momento se siguen de:

$latex frac{partial}{partial x^nu} T^{i nu} = 0$

que es:

$latex frac{d}{dt} q = – frac{1}{N} nabla P$

y en forma SPH queda:

$latex 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 $latex nu_b$ es el número de bariones asociados a la partícula $latex b$.

Y la de la energía se sigue de:

$latex frac{partial}{partial x^j} T^{0j} = 0$

que es:

$latex frac{d}{dt} epsilon = – frac{1}{N} nabla cdot (Pv)$

y en foma SPH queda:

$latex 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} $

Tags: , , , , , ,

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 $latex A(r)$ de la siguiente manera:

$latex A_I(r) = int A(r’) W(r-r’,h) dr’$

donde $latex r$ es el vector posición, $latex W$ es una función de pesos o kernel (con las propiedades ya comentadas en otros posts) y $latex 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:

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

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

Tags: , ,

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, $latex frac{d}{dt}$, se relaciona con la derivada Euleriana respecto al tiempo, $latex frac{partial}{partial t}$ de la siguiente manera:

$latex 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: $latex frac{d}{dt} rho = – rho nabla cdot vec{v}$
  2. Ecuacion del momento: $latex frac{d}{dt} vec{v} = -frac{nabla P}{rho} + vec{f}$
  3. Ecuación de la energía: $latex 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: $latex P = (gamma -1) cdot rho cdot epsilon$ (ecuación del gas ideal)

Tags: , , , , , , , , , , ,

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

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

tienen parte real $latex 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.

 

Tags: , , , ,

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.

Tags: ,

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

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

Tags: , , , ,

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:

[sourcecode languaje=”cpp”]

#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();

};
[/sourcecode]

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

[sourcecode languaje=”cpp”]
#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 …;
}

};
[/sourcecode]

Una de las funcionalidades soportadas a partir de Fortran 2003 es la orientación a objetos. De esta manera, tenemos herencia, abstracción y polimorfismo en el lenguaje para computación de altas prestaciones por excelencia.

Sin embargo, cuando hablamos de Fortran 2003 nos referimos al estandar correspondiente, cuya implementación corresponde a los diferentes compiladores.

En particular, el compilador gfortran de GNU, ¿qué funcionalidades incorpora del estardar 2003? ¿qué es lo que incorpora de la Orientación a Objetos? Y de la misma manera, ¿qué soportan el resto de compiladores? En este enlace tenemos información al respecto.

Tags: , , , , , ,

En el artíclo [Gingold & Monaghan, 1977] se introducen

La función kernel, $latex W(vec{r}-vec{r’},h)$, es una función que permite interpolar los valores de cualquier propiedad del fluido en función del valor de las partículas del entorno. Su papel es similar al de los diferentes esquemas de diferencias en el ámbito de las Diferencias Finitas o las funciones de forma en los Elementos Finitos.

Existen diferentes funciones kernel: Gaussiana, cuadrática, spline cúbico, quíntica, etc.

La función kernel debe cumplir:

  1. Positiva: $latex W(r-r’,h) geq 0$ dentro del dominio.
  2. Soporte compacto: $latex W(r-r’,h) = 0$ fuera del dominio.
  3. Normalizada: $latex int W(r-r’,h) dr’ = 1$.
  4. Comportamiento de función delta: $latex lim_{h rightarrow 0} W(r-r’,h) dr’ = delta(r-r’)$.
  5. Monotona decreciente.

Para derivar, tomamos la derivada analítica de la suma aproximada. De esta manera, como la derivada de la función kernel es conocida, no necesitamos diferencias finitas y el conjunto de ecuaciones PDE pasa a ser ODE.

$latex nabla f(vec{r}) = sum_b frac{m_b}{rho_b} f_b W(vec{r}-vec{r’},h)$

Tags: , , , ,

En el artículo  [Gingold & Monaghan, 1977] se presenta por primera vez el método Smoothed Particle Hydrodynamics. Los autores, originalmente, buscaban un método que permitiera tratar problemas en astrofísica asimétricos (sin simetría esférica, sin simetría axial, etc.) . En estos casos, los métodos de diferencias finitas no se adaptan bien, pues requieren elevar el número de puntos en la malla para seguir con la precisión deseada su evolución, lo cual complica enormemente la evolución de las integrales múltiples.

Lo que pensaron es utilizar la descripción Lagrangiana del flujo del fluido que centra su atención en los elementos del fluido. En la discretización, estos elementos se mueven siguiendo las leyes de Newton con fuerzas debidas a los gradientes de presión y a otras fuerzas de cuerpo como la gravedad, rotación o magnéticas. ¿Qué método utilizar para determinar las fuerzas que actuan en un momento determinado sobre un elemento de fluido?

Para empezar, para elementos de fluido de igual masa, el número de elementos por unidad de volumen debe ser proporcional a la densidad. Además, sin ningún tipo especial de simetría, la posicion de los elementos será aleatoria de acuerdo con la densidad. Para recuperar la densidad de la distribución conocida de los elementos es equivalente a recuperar una distribución de probabilidad a partir de una muestra. Existen dos métodos para conseguir esto que funcionan bien con problemas de fluidos: el smoothing kernel method y  la técnica del spline delta. Ambos métodos se pueden pensar como la aproximación de una integral por el procedimiento de Monte Carlo.

Tags: , , , , , , , ,

En relatividad especial, o relatividad restringida, podemos generalizar las ecuaciones de Euler clásicas que determinan la evolución de los fluidos.

Para ello, definimos el tensor de energia-impulso de la siguiente manera:

$latex T^{alpha beta} = (rho + frac{p}{c^2}) u^alpha u^beta + p eta^{alpha beta}$

De esta manera:

$latex frac{partial}{partial x^beta} T^{alpha beta} = 0$

generaliza tanto la ecuación de conservación de la masa como las ecuaciones de Euler.

Tags: , , , ,

Las estrellas con una masa mayor a 8 masas solares finalizan su evolución hidroestática con una supernova. El remanente de este evento puede ser tanto un agujero negro, en el caso de las estrellas mas masivas , o una estrella compacta, comunmente llamada estrella de neutrones, para el resto de estrellas.

Existe la posibilidad de que estos objetos sean en realidad estrellas de quarks.

Tags: , , , ,

En arquitectura de computadores, la taxonomía de Flynn clasifica a las mismas en función del número de instrucciones que se pueden ejecutar de forma concurrente y el número de datos por instrucción.

Una de las posibilidades es la arquitectura  SIMD (Single Instruction Multiple Data), que explota el paralelismo a nivel de datos, puesto que cada instruccion se puede ejecutar simultaneamente sobre un conjunto amplio de datos.

Las GPU (Graphics Processing Unit) son procesadores desarrollados para el procesamiento gráico que implementan una arquitectura vectorial. Si queremos utilizar de esta potencia para nuestros programas paralelos, disponemos de  CUDA y OpenCL.

CUDA (Compute Unified Device Architecturees) es el sofware que ofrece NVIDIA para explotar la potencia de las GPUs en las aplicaciones programadas por el usuario.

OpenCL (Open Computing Language) es la alternativa libre a CUDA. No solo pretende trabajar con GPUs sino también con CPUs. La arquitectura OpenCL se muestra en el siguiente gráfico:

Tags: , , , , , ,

Las ecuaciones de Euler gobiernan la dinámica de los fluidos compresibles, como gases o líquidos a alta presión, cuando consideramos despreciables las fuerzas de cuerpo, las tensiones viscosas y los flujos de calor. Forman un sistema de PDE no lineal hiperbólico.

En el caso clásico, deducidas por Leonhard Euler, las leyes de conservación son las siguientes:

  1. Conservación de la masa: $latex rho_t + nabla cdot (rho vec{v}) = 0$
  2. Conservación del momento: $latex (rho vec{v})_t + nabla cdot (rho vec{v} otimes vec{v} + pI) = 0 $
  3. Conservación de la energia: $latex E_t + nabla cdot [(E + p)] vec{v} = 0$

donde $latex rho(x,y,z,t)$ es la densidad de masa, $latex vec{v}= (v^1,v^2,v^3)$ es el vector velocidad con $latex v^i(x,y,z,t)$, $latex p(x,y,z,t)$ es la presión,

De hecho, podemos escribir el sistema de forma compacta:

$latex U_t + nabla cdot H = 0$

con el vector columna $latex U = left[ begin{array}{c} rho \ rho v^1 \ rho v^2 \ rho v^3 \ E end{array} right]$ y el tensor $latex H = begin{bmatrix} rho v^1 & rho (v^1)^2 + p & rho v^2 v^1 & rho v^3 v^1 & v^1 (E + p) \ rho v^2 & rho v^1 v^2 & rho (v^2)^2 + p & rho v^3 v^2 & v^2 (E + p) \ rho v^3 & rho v^1 v^3 & rho v^2 v^3 & rho (v^3)^2 + p & v^3 (E + p) end{bmatrix}$

La derivación de estas leyes de conservación esta basada en la relación entre las integrales en volumenes de control y sus fronteras y utilizando el teorema de Gauss. En la forma integral de las ecuaciones no necesitamos la hipótesis de diferenciabilidad. Para la conservación de la masa asumimos que en un volumen $latex V$ la masa ni se crea ni se destruye, por lo que la variación de fluido en su interior esta relacionada con la cantidad del mismo que atraviesa su frontera $latex partial V$

Tags: , , , ,

La hidrodinámica (HD) es la parte de la física que estudia la dinámica de los fluidos tanto incompresibles, los líquidos, como compresibles, los gases o los líquidos a alta presión (de hecho, todos los fluidos son compresibles, siendo la incompresibilidad una aproximación para simplificar las ecuaciones que describen su dinámica).

La magnetohidrodinámica (MHD) estudia la dinámica de fluidos conductores de electricidad en presencia de campos electromagnéticos. El conjunto de ecuaciones que describen la MHD son una combinación de las ecuaciones de Navier Stokes de la dinámica de fluidos y las ecuaciones de Maxwell del electromagnetismo que deben ser resueltas simultaneamente.

Cuando tenemos flujos a velocidades cercanas a la velocidad de la luz entonces hablamos de hidrodinámica en relatividad especial (SRHD) y magnetohidrodinámica en relatividad especial (SRMHD).

Finalmente, cuando el fluido está en presencia de fuertes campos gravitatorios, como por ejemplo en presencia de objetos compactos, hablamos de hidrodinámica y magnetohidrodinámica en relatividad general (GRHD y GRMHD).

Tags: , , , , ,

FireStats icon Powered by FireStats