You are currently browsing the tag archive for the ‘Gauss-Seidel’ tag.

Llego ya un tiempo implementando los diferentes esquemas multigrid en C++ y, tras un periodo de larga oscuridad, por fin aparece un poco de luz: ayer , por primera vez, pareció funcionar todo (en una dimensión y con fronteras Dirichlet, claro… Aunque todo lo demás ya está también implementado).

En el primer nivel necesitamos un esquema de relajación. Esta implementado un weighted Jacobi, por su, a mi humilde entender, natural paralelización en CUDA (y que cuando w=1 tenemos un Jacobi), un SOR (sobrerelajación, que con w=1 deriva en un Gauss-Seidel), que suele ser la opción de siempre y un weighted red-black Gauss-Seidel, de cara a paralelizaciones tradicionales (OpenMP, MPI).

En el segundo nivel tenemos los algoritmos de corrección con esquemas V-Cycle, en iterativo, y \mu-Cycle, en recursivo.

Finalmente, tenemos el FMG en el último nivel, tambien en iterativo y en recursivo. En total, 12 combinaciones posibles (en realidad, doce para una dimensión, doce para dos dimensiones y doce para tres dimensiones)

A ver que tal se dan hoy las ejecuciones en dos dimensiones.

Anuncios

Ya comentamos en el post anterior que formato tienen los sistemas Au=f que surgen a partir de la discretización de un problema de Poisson y que tendremos que resolver. Esta resolución la podriamos hacer de manera directa pero son mucho mas eficaces los métodos iterativos.

Se puede demostrar que la manera mas sencilla para obtener los algoritmos de los diferentes métodos iterativos es pensar sencillamente en la aproximación por diferencias finitas, despejando la variable central i calcularla en función de los vecinos del paso de tiempo anterior para Jacobi, introduciendo un peso para ponderar el paso anterior de la misma variable con el calculado por Jacobi para relajación y, finalmente, utilizar no las varialbles del paso anterior sino tambien las ya calculadas del paso actual para Gauss-Seidel.

Además, el hecho de trabajar con n dimensiones solo significa colocar nuestra expresión en n bucles anidados. Así de sencillo…

En el caso de 1D, tendriamos:

for (int i=1; i<nx-1; i++) {

    //Jacobi
    vj[t+1][i] = 0.5 * ( vj[t][i-1] + vj[t][i+1] + h*h*f[i] );

    //weighted Jacobi
    vwj[t+1][i] = vwj[t][i] + w * ( vwj[t][i] - vj[t+1][i] );

    //Gauss-Seidel
    vgs[i] = 0.5 * ( vgs[i-1] + vgs[i+1] + h*h*f[i] );

    //weighted Gauss-Seidel
    vwgs[i] = vwgs[i] + w * ( vwgs[i] - vgs[i] );

}
julio 2018
L M X J V S D
« Ago    
 1
2345678
9101112131415
16171819202122
23242526272829
3031  
Anuncios