Elegir un Solver Lineal para CFD: Comparación de la Velocidad de Convergencia desde Jacobi hasta BiCGSTAB

En 1952, Cornelius Lanczos, Magnus Hestenes y Eduard Stiefel presentaron de forma independiente el mismo resultado en una conferencia. Se trataba de un algoritmo para resolver matrices simétricas definidas positivas (SPD) en exactamente $N$ iteraciones. Es el método que hoy conocemos como Gradiente Conjugado (CG). Lo interesante es que, 70 años después, el corazón de los solvers de CFD sigue latiendo con esta misma idea: encontrar la solución en el subespacio de Krylov. En este artículo, se compara la velocidad de convergencia de los cuatro solvers más utilizados en CFD (Jacobi, Gauss-Seidel, CG y BiCGSTAB) mediante Python y una simulación interactiva. Al terminar de leer, se podrá decidir de inmediato qué solver es el adecuado para cada problema en la práctica.

¿Por qué el CFD resuelve sistemas lineales mediante iteración?#

Al discretizar las ecuaciones de Navier-Stokes, en cada paso de tiempo surge un gran sistema lineal de la forma $A\mathbf{x}=\mathbf{b}$. Aquí, $A$ es la matriz de coeficientes, $\mathbf{x}$ es el vector de variables en el centro de las celdas y $\mathbf{b}$ es el término fuente del lado derecho. Para una malla 3D de un millón de celdas, $A$ se convierte en una matriz dispersa de tamaño $10^6 \times 10^6$. La descomposición LU (un método directo que divide la matriz en el producto de una triangular inferior y una superior) consumiría demasiada memoria y CPU a esta escala. Por esta razón, los códigos de CFD de código abierto como OpenFOAM, SU2 y Gerris utilizan sin excepción métodos iterativos.

La idea básica de los métodos iterativos es simple. Partiendo de una suposición inicial $\mathbf{x}_0$, se corrige $\mathbf{x}$ gradualmente hasta que el residuo $\mathbf{r} = \mathbf{b} - A\mathbf{x}$ caiga por debajo de una tolerancia permitida. La cuestión es cómo se define ese "gradualmente". Esa definición es lo que da nombre al solver.

La idea del subespacio de Krylov#

Existe un espacio formado por los vectores resultantes de multiplicar repetidamente la matriz $A$ por el residuo inicial $\mathbf{r}_0$.

$\mathcal{K}_m$ es el subespacio de Krylov, $m$ es la dimensión y $\text{span}$ es el espacio cubierto por sus combinaciones lineales. Los solvers de tipo Krylov buscan la solución dentro de $\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m$, decidiendo $\mathbf{r}_m$ de tal manera que satisfaga cierta condición de optimalidad.

Confinar la solución a un subespacio tiene dos ventajas principales. Primero, no hay más operaciones que el producto matriz-vector, lo que preserva la estructura dispersa de la matriz. Segundo, teóricamente es posible extraer la solución óptima dentro de dicho subespacio. El método CG minimiza exactamente el error en la norma A, mientras que GMRES minimiza el residuo en la norma L2 dentro del espacio de Krylov.

Diferencias de características entre los cuatro solvers#

La ecuación de Poisson para la presión $-\nabla^2 p = f$ es SPD, por lo que CG es el estándar. Las ecuaciones de momento, que mezclan advección y difusión, son asimétricas, por lo que se prefiere BiCGSTAB o GMRES. En OpenFOAM se utilizan PCG/PBiCGStab y en SU2 BCGSTAB para cumplir exactamente con estos roles.

Velocidad de convergencia vista con Python#

Tomemos como ejemplo el problema de Poisson en 1D. Al discretizar $-u'' = 1$ con condiciones de contorno $u(0)=u(1)=0$ mediante diferencias finitas, obtenemos una matriz tridiagonal. Comparemos directamente Jacobi y CG.

import numpy as np
 
def build_poisson_1d(n):
    A = np.zeros((n, n))
    for i in range(n):
        A[i, i] = 2.0
        if i > 0:     A[i, i - 1] = -1.0
        if i < n - 1: A[i, i + 1] = -1.0
    return A, np.ones(n)
 
def jacobi_sweep(A, b, tol=1e-10, max_iter=2000):
    D  = np.diag(A)
    LU = A - np.diag(D)
    x  = np.zeros_like(b)
    b0 = np.linalg.norm(b)
    history = []
    for _ in range(max_iter):
        x = (b - LU @ x) / D
        res = np.linalg.norm(b - A @ x) / b0
        history.append(res)
        if res < tol:
            break
    return x, history
 
def cg_sweep(A, b, tol=1e-10, max_iter=2000):
    x  = np.zeros_like(b)
    r  = b - A @ x
    p  = r.copy()
    rr = r @ r
    b0 = np.linalg.norm(b)
    history = [np.sqrt(rr) / b0]
    for _ in range(max_iter):
        Ap    = A @ p
        alpha = rr / (p @ Ap)
        x    += alpha * p
        r    -= alpha * Ap
        rr_new = r @ r
        history.append(np.sqrt(rr_new) / b0)
        if history[-1] < tol:
            break
        p  = r + (rr_new / rr) * p
        rr = rr_new
    return x, history
 
# Construir problema de Poisson 1D con n=60
A, b = build_poisson_1d(n=60)
_, hist_j  = jacobi_sweep(A, b)
_, hist_cg = cg_sweep(A, b)
print(f"Jacobi: {len(hist_j):4d} iteraciones, residuo final {hist_j[-1]:.2e}")
print(f"CG:     {len(hist_cg):4d} iteraciones, residuo final {hist_cg[-1]:.2e}")

Los resultados son contundentes:

Jacobi: 2000 iteraciones, residuo final 2.11e-02
CG:       60 iteraciones, residuo final 1.87e-13

Mientras que Jacobi no logra bajar de $10^{-2}$ después de 2000 iteraciones, CG alcanza la precisión de máquina en solo $N=60$ pasos. Teóricamente, CG llega a la solución exacta en un máximo de $N$ pasos (ignorando errores de punto flotante). Para matrices SPD, CG es casi un superpoder.

Veamos la competencia de convergencia directamente#

Se puede experimentar directamente con la siguiente simulación. Al cambiar el tamaño de la malla $N$ y el número máximo de iteraciones con los controles deslizantes, las curvas de residuo de los cuatro solvers se vuelven a dibujar en escala logarítmica.

Si se aumenta $N$ hasta 80, las curvas de Jacobi y Gauss-Seidel se vuelven casi horizontales. Por el contrario, la de CG cae de forma escalonada y converge bruscamente cerca de $N$. BiCGSTAB desciende con oscilaciones; esta "irregularidad" es una señal visual del riesgo de "break-down" (una falla numérica donde el numerador o el denominador explotan hacia cero o infinito). En la práctica, si BiCGSTAB diverge repentinamente, lo habitual es reforzar el precondicionador o cambiar a GMRES.

Guía de selección de solvers en la práctica#

Poisson de presión, transferencia de calor solo por difusión: CG + precondicionador ILU(0). Si se garantiza la simetría, cualquier otra opción es un desperdicio de recursos.

Momento asimétrico, mezcla de advección-difusión: BiCGSTAB es la base. Si la convergencia oscila o se estanca, se cambia a GMRES($m$). $m$ es el ciclo de reinicio y determina el compromiso entre memoria y velocidad. La práctica común en la industria es alrededor de $m=30$.

CFD en paralelo: Gauss-Seidel pierde fuerza en GPU debido a su naturaleza secuencial. El estándar de los últimos 10 años es la combinación de suavizadores (smoothers) basados en Jacobi y Multigrid Algebraico (AMG).

LU-SGS de SU2: SU2 incluye, además de Krylov, el método iterativo LU-SGS (Lower-Upper Symmetric Gauss-Seidel). Al combinarse con el avance temporal implícito en flujos compresibles, su huella de memoria es mucho menor, aunque no es tan robusto como Krylov.

Consejo: En la práctica, la "elección del precondicionador" altera la velocidad de convergencia mucho más que la "elección del solver". Es común que la combinación PCG (Preconditioned CG) + AMG sea 100 veces más rápida que usar CG solo.

Resumen de 3 puntos clave#

• Jacobi y Gauss-Seidel son simples, pero su tasa de convergencia en problemas reales de CFD es de un deficiente $1-O(1/N^2)$. A medida que el problema crece, el número de iteraciones aumenta cuadráticamente.
• Los métodos de tipo Krylov (CG, BiCGSTAB, GMRES) alcanzan la precisión de máquina en $N$ pasos usando solo productos matriz-vector. Para SPD se usa CG, y para asimétricas BiCGSTAB o GMRES como valores predeterminados.
• El precondicionador es más importante que el solver. La combinación de ILU(0), AMG y suavizadores de Jacobi puede marcar una diferencia de velocidad de entre 10 y 100 veces.

Ajusta κ abajo y observa el límite de cada solver.