julio 2013

You are currently browsing the monthly archive for julio 2013.

El software Chombo, del Berkeley Lab, combina los métodos en diferencias finitas con las mallas adaptativas (AMR) para resolver, entre otras, PDEs elípticas.

Las siguientes imágenes, en 2D y 3D, se han obtenido a partir de su AMRPoisson:

schSolSouschSol

Tags: , , ,

Desde mi modesto blog, una pequeña ayuda a esta página entusiasta de la recuperación del uso del futuro de subjuntivo. Os adjunto el enlace y espero le dieren clic y pulsaren like 🙂

https://www.facebook.com/pages/Amigos-del-futuro-subjuntivo/448669705230946

Tags: ,

Sorpresa máxima al abrir hoy el Reader de WordPress: ¡Tao hablando de la hipótesis de Riemann en su blog!. Deja claro al principio del post que simplemente va a poner junto todo lo que se conoce al respecto, sin aportar nada nuevo. ¡Pero es Tao y es la hipótesis de Riemann! ¿Va a quedar sólo en eso? Expectación máxima. ¡Dios nos coja confesados! 😀

Tags: ,

Definimos a los kernels como funciones del tipo:

$latex W_{ab}=W(boldsymbol{r}_a – boldsymbol{r}_b,h)$,

donde $latex a$ es la partícula en la que está centrada la función y $latex b$ es una partícula dentro del soporte compacto de la función kernel, controlado éste último por $latex h$, la smoothing length (longitud de suavizado).

En este post básicamente pretendo aclarar lo que significa $latex nabla_a W_{ab}$ cuando, por ejemplo, tenemos definido el kernel como:

$latex W(q) = alpha_D exp (-q^2)$ con $latex 0 leq q leq 2$.

Para empezar, $latex alpha_D$ es una constante de dimensionalidad, por lo que la fórmula está escrita de manera compacta y sirve para cualquier dimensión. Además, tenemos que $latex q = frac{r}{h}$, siendo $latex r$ la distancia ente las partículas, por lo que:

$latex r = |boldsymbol{r}_a – boldsymbol{r}_b| =^{(3D)} sqrt{(x_a-x_b)^2 + (y_a-y_b)^2 + (z_a-z_b)^2}$.

Si fijamos la posición de la partícula $latex a$, la función que nos da la distancia de esta a cualquier punto dentro del soporte compacto es:

$latex r_a (boldsymbol{r}) = |boldsymbol{r}_a-boldsymbol{r}| =^{(3D)} sqrt{(x_a-x)^2 + (y_a-y)^2 + (z_a-z)^2}$,

siendo $latex q_a$ lo mismo añadiendo el factor $latex h$.

Por lo tanto, en este caso tenemos, en tres dimensiones y donde $latex b$ es una partícula en una posición arbitraria $latex (x,y,z)$:

$latex nabla_a W_{ab}(q) =^{(3D)} (partial_x W_{ab}(q_a), partial_y W_{ab}(q_a), partial_z W_{ab}(q_a)) =$

$latex = alpha_D exp(-q^2) (-2q) (partial_x (q_a), partial_y (q_a), partial_z (q_a)) = alpha_D exp(-q^2) (-2q) nabla_a q_a$

donde:

$latex nabla_a q_a = frac{-1}{h r_a} (x_a-x,y_a-y,z_a-z)$.

De la misma manera, si tenemos:

$latex W(q) = alpha_D begin{cases} 1-frac{3}{2} q^2 + frac{3}{4} q^3, 0 leq q < 1\ frac{1}{4} (2-q)^3, 1 leq q < 2 \ 0, q geq 2 end{cases}$

entonces:

$latex nabla_a W_{ab}(q) = alpha_D begin{cases} (-3q + frac{9}{4}q^2) nabla_a q_a, 0 leq q < 1 \ -frac{3}{4} (2-q)^2 nabla_a q_a, 1 leq q < 2 \ 0, q geq 2 end{cases}$

Así pues:

$latex nabla W(q) = frac{d}{dq} W(q) nabla q$.

 

Tags: , , , ,

A la hora de resolver las diferentes ecuaciones elípticas CFC tenemos dos posibilidades para fijar las condiciones en la frontera, cada una con sus mas y sus menos.

La primera consiste en hacer un desarrollo multipolar de los terminos fuente en armónicos esféricos, de manera que cuantos mas términos consideremos mas cerca podremos colocar la frontera.

La segunda consiste en compactificar el dominio, lo que nos permite reducir todo el universo a un cubo unidad y considerar Minkowski en su frontera, puesto que ésta corresponde a infinito.

Tags: , ,

En la discretización que hicimos teníamos dos sistemas acoplados, uno para las $latex X^i$ y otro para las $latex beta^i$. Procedemos ahora a desacoplarlos.

Para empezar, tomamos la divergencia (plana) del sistema:

$latex Delta X^i = 8 pi f^{ij} S^*_j – frac{1}{3}mathcal{D}^i mathcal{D}_j X^j$

y, teniendo en cuenta que $latex mathcal{D}$ conmuta con $latex Delta$ (métrica plana), tenemos:

$latex Delta (mathcal{D}_i X^i) = 8 pi mathcal{D}^j S^*_j – frac{1}{3} Delta (mathcal{D}_j X^j)$,

por lo que:

$latex Delta (mathcal{D}_i X^i) = frac{3}{4} 8 pi mathcal{D}^j S^*_j$.

De esta manera, si definimos $latex Theta_X := mathcal{D}_i X^i$, nos queda:

$latex Delta Theta_X = frac{3}{4} 8 pi mathcal{D}^j S^*_j = 6 pi (partial_x S^*_x + partial_y S^*_y +partial_z S^*_z )$,

que discretizado queda:

$latex frac{(Theta_X)_{i-1,j,k}-2(Theta_X)_{i,j,k}+(Theta_X)_{i+1,j,k}}{h_x^2} + $

$latex frac{(Theta_X)_{i,j-1,k}-2(Theta_X)_{i,j,k}+(Theta_X)_{i,j+1,k}}{h_y^2} + $

$latex frac{(Theta_X)_{i,j,k-1}-2(Theta_X)_{i,j,k}+(Theta_X)_{i,j,k+1}}{h_z^2} = $

$latex = 6 pi (partial_x S^*_x + partial_y S^*_y +partial_z S^*_z )_{i,j,k} $,

donde inicialmente:

$latex (S^*_a)_{i,j,k} = (psi^6)_{i,j,k}rho_{i,j,k}h_{i,j,k}w^2_{i,j,k}(v_a)_{i,j,k}$,

$latex (partial_x S^*_x + partial_y S^*_y +partial_z S^*_z )_{i,j,k} = $

$latex frac{(S^*_x)_{i+1,j,k}-(S^*_x)_{i-1,j,k}}{2h_x} + frac{(S^*_x)_{i,j+1,k}-(S^*_x)_{i,j-1,k}}{2h_y} + frac{(S^*_x)_{i,j,k+1}-(S^*_x)_{i,j,k-1}}{2h_z}$

y que es lineal.

El primer sistema acoplado de ecuaciones quedaría ahora:

$latex partial_{xx} X^x + partial_{yy} X^x + partial_{zz} X^x = 8 pi S^*_x – frac{1}{3} partial_x Theta_X approx $

$latex approx frac{X^x_{i-1,j,k}-2X^x_{i,j,k}+X^x_{i+1,j,k}}{h_x^2} + frac{X^x_{i,j-1,k}-2X^x_{i,j,k}+X^x_{i,j+1,k}}{h_y^2} + frac{X^x_{i,j,k-1}-2X^x_{i,j,k}+X^x_{i,j,k+1}}{h_z^2} = $

$latex = 8 pi (S^*_x)_{i,j,k} – frac{1}{3} (partial_x Theta_X)_{i,j,k}$,

¡que vuelve a ser lineal!

Continuamos con:

$latex partial_{xx} X^y + partial_{yy} X^y + partial_{zz} X^y = 8 pi S^*_y – frac{1}{3} partial_y Theta_X approx $

$latex approx frac{X^y_{i-1,j,k}-2X^y_{i,j,k}+X^y_{i+1,j,k}}{h_x^2} + frac{X^y_{i,j-1,k}-2X^y_{i,j,k}+X^y_{i,j+1,k}}{h_y^2} + frac{X^y_{i,j,k-1}-2X^y_{i,j,k}+X^y_{i,j,k+1}}{h_z^2} = $

$latex = 8 pi (S^*_y)_{i,j,k} – frac{1}{3} (partial_y Theta_X)_{i,j,k}$

y, finalmente:

$latex partial_{xx} X^z + partial_{yy} X^z + partial_{zz} X^z = 8 pi S^*_z – frac{1}{3} partial_z Theta_X approx $

$latex approx frac{X^z_{i-1,j,k}-2X^z_{i,j,k}+X^z_{i+1,j,k}}{h_x^2} + frac{X^z_{i,j-1,k}-2X^z_{i,j,k}+X^z_{i,j+1,k}}{h_y^2} + frac{X^z_{i,j,k-1}-2X^z_{i,j,k}+X^z_{i,j,k+1}}{h_z^2} = $

$latex = 8 pi (S^*_z)_{i,j,k} – frac{1}{3} (partial_z Theta_X)_{i,j,k}$,

donde calculamos al principio:

$latex (partial_x Theta_X)_{i,j,k} = frac{(Theta_X)_{i+1,j,k}-(Theta_X)_{i-1,j,k}}{2h_x}$

$latex (partial_y Theta_X)_{i,j,k} = frac{(Theta_X)_{i,j+1,k}-(Theta_X)_{i,j-1,k}}{2h_y}$

$latex (partial_z Theta_X)_{i,j,k} = frac{(Theta_X)_{i,j,k+1} – (Theta_X)_{i,j,k-1}}{2h_z}$

A continuación, discretizamos las siguientes ecuaciones:

$latex hat{A}^{xx} = 2 partial_x X^x – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z) approx$

$latex approx frac{2}{3}frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} -frac{1}{3} frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} – frac{1}{3} frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = hat{A}^{xx}_{i,j,k}$,

$latex hat{A}^{xy} = hat{A}^{yx}= partial_x X^y + partial_y X^x approx $

$latex approx frac{X^y_{i+1,j,k}-X^y_{i-1,j,k}}{2h_x} + frac{X^x_{i,j+1,k}-X^x_{i,j-1,k}}{2h_y} = hat{A}^{xy}_{i,j,k} = hat{A}^{yx}_{i,j,k}$,

$latex hat{A}^{xz} = hat{A}^{zx} = partial_x X^z + partial_z X^x approx $

$latex approx frac{X^z_{i+1,j,k}-X^z_{i-1,j,k}}{2h_x} + frac{X^x_{i,j,k+1}-X^x_{i,j,k-1}}{2h_z} = hat{A}^{xz}_{i,j,k} = hat{A}^{zx}_{i,j,k}$,

$latex hat{A}^{yy} = 2 partial_y X^y – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z) approx $

$latex approx -frac{1}{3}frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} +frac{2}{3} frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} – frac{1}{3} frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = hat{A}^{yy}_{i,j,k}$,

$latex hat{A}^{yz} = hat{A}^{zy} = partial_y X^z + partial_z X^y approx $

$latex approx frac{X^z_{i,j+1,k}-X^z_{i,j-1,k}}{2h_y} + frac{X^y_{i,j,k+1}-X^y_{i,j,k-1}}{2h_z} = hat{A}^{yz}_{i,j,k} = hat{A}^{zy}_{i,j,k}$,

$latex hat{A}^{zz} = 2 partial_z X^z – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z) approx $

$latex approx -frac{1}{3}frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} -frac{1}{3} frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} + frac{2}{3} frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = hat{A}^{zz}_{i,j,k}$.

Por tanto, la siguiente ecuación:

$latex Delta psi = -2 pi psi^{-1} E^* – psi^{-7} frac{(hat{A}^{xx})^2+(hat{A}^{yy})^2+(hat{A}^{zz})^2+2(hat{A}^{xy})^2+2(hat{A}^{xz})^2+2(hat{A}^{yz})^2}{8}$

queda:

$latex approx frac{psi_{i-1,j,k}-2psi_{i,j,k}+psi_{i+1,j,k}}{h_x^2} + frac{psi_{i,j-1,k}-2psi_{i,j,k}+psi_{i,j+1,k}}{h_y^2} + frac{psi_{i,j,k-1}-2psi_{i,j,k}+psi_{i,j,k+1}}{h_z^2} = $

$latex =-2 pi psi^{-1}_{i,j,k} E^*_{i,j,k} – $

$latex – frac{psi^{-7}_{i,j,k}}{8} ( (hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2 ) $,

con:

$latex partial_{psi_{i,j,k}} F(psi_{i,j,k}) = -2 ( frac{1}{h_x^2} + frac{1}{h_y^2} + frac{1}{h_z^2} ) -2 pi psi_{i,j,k}^{-2} E^*_{i,j,k} – $

$latex – frac{7}{8} psi^{-8}_{i,j,k} ( (hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2 )$,

donde:

$latex E^*_{i,j,k} = psi^{6}_{i,j,k} (D_{i,j,k}+tau_{i,j,k})$

y la ecuación:

$latex Delta (alphapsi) = (alpha psi) (2 pi psi^{-2} (E^*+2S^*) + $

$latex + frac{7}{8} psi^{-8} ((hat{A}^{xx})^2+(hat{A}^{yy})^2+(hat{A}^{zz})^2+2(hat{A}^{xy})^2+2(hat{A}^{xz})^2+2(hat{A}^{yz})^2) )$

como:

$latex approx frac{(alphapsi)_{i-1,j,k} – 2(alphapsi)_{i,j,k}+(alphapsi)_{i+1,j,k}}{h_x^2} + $

$latex + frac{(alphapsi)_{i,j-1,k}-2(alphapsi)_{i,j,k}+(alphapsi)_{i,j+1,k}}{h_y^2} + $

$latex + frac{(alphapsi)_{i,j,k-1}-2(alphapsi)_{i,j,k}+(alphapsi)_{i,j,k+1}}{h_z^2} = $

$latex = (alpha psi)_{i,j,k} (2 pi psi^{-2}_{i,j,k} (E^*_{i,j,k}+2S^*_{i,j,k}) + $

$latex + frac{7}{8} psi^{-8}_{i,j,k} ((hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2) )$,

donde:

$latex partial_{(alpha psi)_{i,j,k}} F((alpha psi)_{i,j,k}) = -2 ( frac{1}{h_x^2} + frac{1}{h_y^2} + frac{1}{h_z^2} ) – 2 pi psi^{-2}_{i,j,k} (E^*_{i,j,k}+2S^*_{i,j,k}) + $

$latex – frac{7}{8} psi^{-8}_{i,j,k} ((hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2) )$

con:

$latex S^*_{i,j,k} = psi^6_{i,j,k}(rho_{i,j,k}h_{i,j,k}(w^2_{i,j,k}-1) + 3 p_{i,j,k})$.

Finalmente, tenemos el otro sistema acoplado:

$latex Delta beta^i = mathcal{D}_j(2 alpha psi^{-6} hat{A}^{ij}) – frac{1}{3} mathcal{D}^i(mathcal{D}_j beta^j)$,

con el que procedemos de igual manera que con las $latex X^i$:

$latex Delta(mathcal{D}_i beta^i) = mathcal{D}_i (mathcal{D}_j (2 alpha psi^{-6} hat{A}^{ij})) – frac{1}{3} Delta (mathcal{D}_i beta^i)$,

de manera que:

$latex Delta Theta_beta = frac{3}{4} mathcal{D}^i (mathcal{D}_j (2 alpha psi^{-6} hat{A}^{ij})) =$

$latex frac{3}{2}(partial_{xx}(alpha psi^{-6} hat{A}^{xx}) + partial_{yy}(alpha psi^{-6} hat{A}^{yy}) + partial_{zz}(alpha psi^{-6} hat{A}^{zz})$,

con:

$latex Theta_beta := mathcal{D}_i beta^i$,

que discretizada queda:

$latex frac{(Theta_beta)_{i-1,j,k}-2(Theta_beta)_{i,j,k}+(Theta_beta)_{i+1,j,k}}{h_x^2} + $

$latex frac{(Theta_beta)_{i,j-1,k}-2(Theta_beta)_{i,j,k}+(Theta_beta)_{i,j+1,k}}{h_y^2} + $

$latex frac{(Theta_beta)_{i,j,k-1}-2(Theta_beta)_{i,j,k}+(Theta_beta)_{i,j,k+1}}{h_z^2} = $

$latex frac{3}{2}((partial_{xx}(alpha psi^{-6} hat{A}^{xx}))_{i,j,k} + (partial_{yy}(alpha psi^{-6} hat{A}^{yy}))_{i,j,k} + (partial_{zz}(alpha psi^{-6} hat{A}^{zz})_{i,j,k})$,

De esta manera, tenemos:

$latex Delta beta^x = partial_x (2 alpha psi^{-6} hat{A}^{xx}) + partial_y (2 alpha psi^{-6} hat{A}^{xy}) + partial_z (2 alpha psi^{-6} hat{A}^{xz}) – frac{1}{2} partial_x Theta_beta approx $

$latex approx frac{beta^x_{i-1,j,k}-2beta^x_{i,j,k}+beta^x_{i+1,j,k}}{h_x^2} + frac{beta^x_{i,j-1,k}-2beta^x_{i,j,k}+beta^x_{i,j+1,k}}{h_y^2} + frac{beta^x_{i,j,k-1}-2beta^x_{i,j,k}+beta^x_{i,j,k+1}}{h_z^2} = $

$latex = (partial_x (2 alpha psi^{-6} hat{A}^{xx}))_{i,j,k} + (partial_y (2 alpha psi^{-6} hat{A}^{xy}))_{i,j,k} + (partial_z (2 alpha psi^{-6} hat{A}^{xz}) )_{i,j,k} – $

$latex – frac{1}{3} (partial_x Theta_beta)_{i,j,k}$.

De la misma manera:

$latex Delta beta^y = partial_x (2 alpha psi^{-6} hat{A}^{yx}) + partial_y (2 alpha psi^{-6} hat{A}^{yy}) + partial_z (2 alpha psi^{-6} hat{A}^{yz}) – frac{1}{3} partial_y Theta_beta approx $

$latex approx frac{beta^y_{i-1,j,k}-2beta^y_{i,j,k}+beta^y_{i+1,j,k}}{h_x^2} + frac{beta^y_{i,j-1,k}-2beta^y_{i,j,k}+beta^y_{i,j+1,k}}{h_y^2} + frac{beta^y_{i,j,k-1}-2beta^y_{i,j,k}+beta^y_{i,j,k+1}}{h_z^2} = $

$latex = (partial_x (2 alpha psi^{-6} hat{A}^{yx}))_{i,j,k} + (partial_y (2 alpha psi^{-6} hat{A}^{yy}))_{i,j,k} + (partial_z (2 alpha psi^{-6} hat{A}^{yz}) )_{i,j,k} – $

$latex – frac{1}{3} (partial_y Theta_beta)_{i,j,k}$.

Y, por último:

$latex Delta beta^z = partial_x (2 alpha psi^{-6} hat{A}^{zx}) + partial_y (2 alpha psi^{-6} hat{A}^{zy}) + partial_z (2 alpha psi^{-6} hat{A}^{zz}) – frac{1}{3} partial_z Theta_beta approx $

$latex approx frac{beta^z_{i-1,j,k}-2beta^z_{i,j,k}+beta^z_{i+1,j,k}}{h_x^2} + frac{beta^z_{i,j-1,k}-2beta^z_{i,j,k}+beta^z_{i,j+1,k}}{h_y^2} + frac{beta^z_{i,j,k-1}-2beta^z_{i,j,k}+beta^z_{i,j,k+1}}{h_z^2} = $

$latex = (partial_x (2 alpha psi^{-6} hat{A}^{zx}))_{i,j,k} + (partial_y (2 alpha psi^{-6} hat{A}^{zy}))_{i,j,k} + (partial_z (2 alpha psi^{-6} hat{A}^{zz}) )_{i,j,k} – $

$latex – frac{1}{3} (partial_z Theta_beta)_{i,j,k}$.

Parece que, del sistema no lineal acoplado inicial, hemos llegado a un sistema de diez ecuaciones desacopladas donde ocho de ellas son lineales y solo dos son no linales. No pinta mal. Ya escribiremos próximamente sobre las condiciones de contorno…

Tags: , , , ,

FireStats icon Powered by FireStats