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

En 1952, Lanczos, Hestenes y Stiefel llegaron al mismo congreso con el mismo algoritmo, escrito por separado. Matrices simétricas definidas positivas (SPD), terminadas en $N$ iteraciones. El nombre que quedó fue gradiente conjugado (CG). Setenta años después, la ecuación de presión de un código CFD sigue buscando la respuesta en ese subespacio.

Al mover el tamaño de malla N de 20 a 80, las cuatro curvas de residuo se dibujan de nuevo. Al subir el tope de iteraciones, cambia el orden en que tocan el suelo.

Cerca de N=40, CG se quiebra hacia abajo y Jacobi queda casi horizontal en el eje logarítmico. Al empujar N hasta 80, la separación crece. Los dientes de BiCGSTAB son un residuo que no muere de forma lisa.

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

Al discretizar Navier–Stokes, cada paso deja un sistema $A\mathbf{x}=\mathbf{b}$. $A$ es la matriz de coeficientes, $\mathbf{x}$ las incógnitas en el centro de celda, $\mathbf{b}$ el término derecho. Con un millón de celdas, $A$ es una matriz dispersa de $10^6 \times 10^6$. La factorización LU (partir en factores triangular inferior y superior) quema memoria y tiempo a la vez en ese tamaño. OpenFOAM, SU2 y Gerris iteran por eso.

Se parte de una estimación $\mathbf{x}_0$ y se corrige $\mathbf{x}$ hasta que el residuo $\mathbf{r} = \mathbf{b} - A\mathbf{x}$ cae bajo la tolerancia. La definición de "poco a poco" es el nombre del solver.

La idea del subespacio de Krylov#

Los productos repetidos de $A$ con el residuo inicial $\mathbf{r}_0$ llenan un espacio.

$\mathcal{K}_m$ es el subespacio de Krylov, $m$ su dimensión, $\text{span}$ el conjunto de combinaciones lineales. La actualización se queda en $\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m$, y el residuo que queda, $\mathbf{r}_m$, tiene que cumplir una condición de optimalidad.

La única operación necesaria es un producto matriz-vector, así que el patrón disperso sobrevive. CG minimiza ahí el error en norma A. GMRES minimiza el residuo en norma 2. Lo segundo ocurre en el subespacio que Arnoldi construye al apilar vectores.

Diferencias de características entre los cuatro solvers#

El Poisson de presión $-\nabla^2 p = f$ es SPD, así que encaja CG. La cantidad de movimiento con advección es no simétrica, así que la elección pasa a BiCGSTAB o GMRES. PCG/PBiCGStab de OpenFOAM y BCGSTAB de SU2 son ese reparto, escrito en el nombre de la opción.

Velocidad de convergencia vista con Python#

El problema modelo es el Poisson 1D, $-u'' = 1$, con $u(0)=u(1)=0$. Las diferencias finitas dejan una matriz tridiagonal. En $N=60$, Jacobi y CG comparten la tolerancia $10^{-10}$.

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}")

La corrida imprime esto.

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

Al subir el control del logaritmo del número de condición, las cuatro curvas se detienen en otro sitio dentro del mismo presupuesto de iteraciones. El tope de iteraciones decide dónde se corta el gráfico.

Con el valor base $\kappa=10^3$ y unas 80 iteraciones, Jacobi (rojo) queda casi horizontal y CG (cian) baja más empinado. Al empujar $\kappa$ hacia $10^6$, Jacobi sigue sin suelo aunque se le den 200 iteraciones. El caso impreso con $N=60$ es un punto de esa figura: Jacobi en 2.11e-2 tras 2000 pasadas, CG en 1.87e-13 tras 60.

El residuo impreso se lee como un factor por paso. De 1 a 2.11e-2 en 2000 pasadas sale una contracción de 0.998 por pasada, la raíz 2000 de 0.0211. Es la escala del radio espectral $1-O(1/N^2)$ de la tabla. Con $N=60$, $1/N^2$ vale 1/3600, cerca de 0.00028, y 0.998 es ese término con una constante pequeña delante.

Si el factor se mantiene, $10^{-10}$ cuesta unas 12.000 iteraciones. Desde 2.11e-2 faltan unas 10.000 pasadas. Al duplicar $N$, $1/N^2$ cae a la cuarta parte y el mismo residuo cuesta cuatro veces más iteraciones.

Las 60 iteraciones de CG igualan a $N$: el techo de como máximo $N$ pasos se gasta. 1.87e-13 es unas $1.1 imes 10^{11}$ veces menor que 2.11e-2. En doble precisión se suma a 1 del orden de $10^{-16}$, y se gastaron dos o tres dígitos. Un millón de celdas a 8 bytes son 8 MB por vector: 32 MB en los cuatro de CG y 56 MB en los siete de BiCGSTAB. Guardar el residuo es barato. La cola de iteraciones no.

La curva con el precondicionador apagado#

Esas curvas son iteraciones desnudas, sin precondicionador. Al encenderlo, una marca del eje horizontal significa otra cosa. Un pico vertical de BiCGSTAB es un break-down: un producto interno del denominador se pegó a cero. Sumar iteraciones no deshace eso. Se cambia el precondicionador, o se pasa a GMRES. Un periodo de reinicio $m$ cerca de 30 es el compromiso habitual entre memoria y ese salto.

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

Poisson de presión, y transferencia de calor que es solo difusión: CG con ILU(0). Si la simetría sigue intacta, otro nombre solo guarda más vectores.

La cantidad de movimiento no simétrica arranca en BiCGSTAB. Cuando el residuo oscila, se cambia a GMRES($m$). Aquí $m$ es el periodo de reinicio.

Gauss-Seidel es secuencial, así que en una GPU se queda sin fuerza. El patrón de la última década es un suavizador de Jacobi más multigrid algebraico (AMG).

SU2 también lleva LU-SGS (Lower-Upper Symmetric Gauss-Seidel), una iteración estacionaria al lado de Krylov. Atado al avance temporal implícito en flujo compresible, la huella de memoria queda más chica. No es tan firme como Krylov.

El precondicionador mueve el número de iteraciones más que el nombre del solver. PCG con AMG suele ser 100 veces más rápido que CG solo. Ese factor es si se corta o no la cola de diez mil pasadas.

Tomar 2.11e-2 por un residual ya convergido#

Entregar al paso de tiempo siguiente un residual de Poisson que sigue en 2.11e-2 tras 2000 pasadas mete ese error de presión en la iteración no lineal que viene. Mientras la matriz siga siendo SPD, el lugar es de CG. En cuanto la advección la vuelve no simétrica, CG deja de ser elegible y el relevo es BiCGSTAB. Elegir el nombre primero y medir el número de condición después es el fallo.

La pregunta que queda es si hace falta apilar otras diez mil iteraciones. Un ciclo V que barre el error en una malla gruesa puede cortar esta cola $1-O(1/N^2)$, o puede volver a llamar a Jacobi como suavizador de la malla fina. Este residual no lo decide.