Donde Jacobi se muere en una malla de capa límite — el precondicionador linelet

La misma ecuación. El mismo número de celdas, la misma iteración, el mismo criterio de convergencia. Bastó aplastar la malla una vez contra la pared para que el número de iteraciones saltara de 18 a 5.000.

La física no cambió en nada. Lo único que cambió es la relación de aspecto de la celda. Y aun así el solver lineal empieza a comportarse como si le hubieran entregado otro problema. Tampoco hay forma de esquivarlo: alcanzar $y^+ \approx 1$ en un cálculo turbulento obliga a aplastar la primera celda junto a la pared.

Ese factor 300 tiene una dirección dentro de la matriz. El precondicionador linelet — una cadena de celdas enhebrada a lo largo de la dirección fuertemente acoplada, el que usa SU2 — lo devuelve. El contador de barridos es la comprobación. Las filas cortas son la tabla. Las filas largas son los coeficientes.

Aplastar una celda convierte la matriz en otra cosa#

Al discretizar el término difusivo por volúmenes finitos sobre una celda rectangular, los coeficientes de cara quedan así.

$\Gamma$ es el coeficiente de difusión y $\Delta x$, $\Delta y$ son las dimensiones de la celda. Ahora el cociente.

$\mathrm{AR}$ es la relación de aspecto de la celda (longitudinal sobre normal a la pared). El cociente de coeficientes va con el cuadrado. Una primera celda con $\mathrm{AR} = 1000$ tiene un acoplamiento norte-sur un millón de veces mayor que el este-oeste. Aplastar la malla diez veces sesga la matriz cien veces.

Al llevar el deslizador AR de 1 a 100, las franjas verticales desaparecen primero. Al activar linelet, la mancha horizontal que queda pierde vida.

Cuando el valor measured bajo el mapa se pega cerca de 0,997, el estancamiento ya está. Con el mismo AR, linelet se lleva esa mancha en uno o dos barridos.

Con la derivada temporal implícita, la diagonal queda

donde $V = \Delta x \Delta y$ es el volumen de la celda y $\Delta t$ el paso de (pseudo-)tiempo. Ese término $V/\Delta t$ es lo único que sostiene la diagonal, y no alcanza el ritmo al que crece $a_N$.

El factor del suavizador puntual lleva AR dentro#

El decaimiento del error de un barrido de Jacobi queda acotado por la suma de los términos fuera de la diagonal sobre la diagonal.

Con $V/\Delta t = 4a_E$ esto vale 0,5 en $\mathrm{AR} = 1$ y 0,9998 en $\mathrm{AR} = 100$. El problema es que la cota se pega a 1 cuando $\mathrm{AR}$ crece y deja de decir nada. Para ver qué sobrevive de verdad hay que ir modo a modo. Para un modo de Fourier $(\theta_x, \theta_y)$ el factor de amortiguamiento es

El modo más lento es el de menor $\theta_y$, es decir, el que varía más suavemente en dirección normal a la pared. Cuando $a_N$ domina la diagonal y $a_P \approx 2a_N$,

donde $N_y$ es el número de capas normales a la pared. Que aquí no aparezca $\mathrm{AR}$ es mala noticia, no buena: seguir aplastando ya no empeora el factor, pero el techo en el que se satura es inservible de por sí. Con $N_y = 48$ vale 0,99786, o sea unos 6.000 barridos para bajar el error un factor $10^{-6}$.

Lo que ocurre se visualiza fácil. Un suavizador puntual borra el error de alta frecuencia normal a la pared en uno o dos barridos. Ahí se atasca. Eliminar el error que varía a lo largo de la pared exige que la información viaje de lado por $a_E$, y como la diagonal está sujeta por $a_N$, cada paso se encoge en $\mathrm{AR}^2$.

Dónde se detiene una iteración que solo mira la diagonal, método por método, está en una comparación de barridos desde Jacobi hasta BiCGSTAB. Allí cambia el solver. Aquí la relación de aspecto de la celda traiciona al mismo solver.

Tres remedios en una sola tabla#

El trato del linelet es explícito: invierte de forma exacta una sola dirección, la fuertemente acoplada. A cambio, el almacenamiento es un único arreglo de tres bandas y el coste cae en el mismo orden de magnitud que Jacobi. Como nunca toca la matriz completa al modo de ILU(0), resulta además mucho menos sensible al particionado MPI.

Cómo se construyen los linelets#

1. En cada celda se mide el cociente de coeficientes de cara $\sigma = a_f / a_{f,\perp}$. En una cara normal a la pared vale $\sigma = \mathrm{AR}^2$.
2. Se marca una cara como fuerte cuando $\sigma \ge \sigma_{\min}$. Un valor por defecto habitual ronda $\sigma_{\min} = 10$.
3. Se arranca en una celda pegada a la pared y se extiende una cadena hacia arriba siguiendo caras fuertes.
4. Se admiten como mucho dos caras fuertes por celda. Si aparece una bifurcación se conserva solo el coeficiente mayor: permitir ramas hace que el resultado deje de ser tridiagonal.
5. La cadena se corta en cuanto alcanza una celda que ya pertenece a otro linelet.
6. Las cadenas de longitud uno se descartan y esas celdas vuelven a Jacobi.

Al renumerar dentro de una cadena, esa submatriz queda tridiagonal. Y tridiagonal significa que el algoritmo de Thomas (eliminación hacia adelante y luego sustitución hacia atrás) la invierte de forma exacta en $O(n)$. Esa frase es todo el precondicionador linelet.

Dónde se corta una cadena queda junto a la tabla de barridos. Los controles son el espesor de la primera capa y la razón de crecimiento.

Contando barridos en Python#

Sobre una malla $64 \times 48$ se deja únicamente error (el término fuente es $f = 0$, así que la solución exacta es $u = 0$), se ejecutan Jacobi puntual y Gauss-Seidel por líneas, y se cuentan los barridos hasta que el error baja $10^{-6}$.

import numpy as np
 
def cell_coefficients(dx, dy, gamma=1.0, diag_factor=4.0):
    """Coeficientes difusivos FVM de una celda rectangular. diag_factor fija V/dt como múltiplo de a_E."""
    a_e = gamma * dy / dx        # caras este/oeste
    a_n = gamma * dx / dy        # caras norte/sur
    d0 = diag_factor * a_e       # término V/dt
    return a_e, a_n, d0
 
def residual_field(u, a_e, a_n, d0):
    r = -d0 * u
    r[1:-1, 1:-1] -= a_e * (2*u[1:-1, 1:-1] - u[:-2, 1:-1] - u[2:, 1:-1])
    r[1:-1, 1:-1] -= a_n * (2*u[1:-1, 1:-1] - u[1:-1, :-2] - u[1:-1, 2:])
    r[0, :] = r[-1, :] = 0.0     # contornos de Dirichlet
    r[:, 0] = r[:, -1] = 0.0
    return r
 
def jacobi_sweep(u, a_e, a_n, d0, omega=1.0):
    r = residual_field(u, a_e, a_n, d0)
    u += omega * r / (d0 + 2*a_e + 2*a_n)
 
def thomas(a, b, c, d):
    """Solución tridiagonal exacta: eliminación hacia adelante y sustitución hacia atrás."""
    n = len(d)
    cp, dp = np.empty(n), np.empty(n)
    cp[0], dp[0] = c[0]/b[0], d[0]/b[0]
    for k in range(1, n):
        m = b[k] - a[k]*cp[k-1]
        cp[k] = c[k]/m
        dp[k] = (d[k] - a[k]*dp[k-1])/m
    x = np.empty(n)
    x[-1] = dp[-1]
    for k in range(n-2, -1, -1):
        x[k] = dp[k] - cp[k]*x[k+1]
    return x
 
def linelet_sweep(u, a_e, a_n, d0):
    """Invierte de forma exacta una cadena normal a la pared completa (Gauss-Seidel por líneas)."""
    nx, ny = u.shape
    m = ny - 2
    dg = d0 + 2*a_e + 2*a_n
    a = np.full(m, -a_n)
    b = np.full(m, dg)
    c = np.full(m, -a_n)
    a[0] = c[-1] = 0.0
    for i in range(1, nx-1):
        rhs = a_e * (u[i-1, 1:-1] + u[i+1, 1:-1])   # el acoplamiento fuera de línea pasa al miembro derecho
        u[i, 1:-1] = thomas(a, b, c, rhs)
 
def count_sweeps(ar, kind, nx=64, ny=48, tol=1e-6, max_sweeps=20000):
    dx, dy = 1.0/nx, 1.0/nx/ar
    a_e, a_n, d0 = cell_coefficients(dx, dy)
    u = np.random.default_rng(7).standard_normal((nx, ny))
    u[0, :] = u[-1, :] = 0.0
    u[:, 0] = u[:, -1] = 0.0
    e0 = prev = np.linalg.norm(u)
    for k in range(1, max_sweeps + 1):
        jacobi_sweep(u, a_e, a_n, d0) if kind == 'jacobi' else linelet_sweep(u, a_e, a_n, d0)
        e = np.linalg.norm(u)
        rho, prev = e/prev, e
        if e/e0 < tol:
            return k, rho
    return max_sweeps, rho
 
print(f"{'AR':>6} {'Jacobi':>8} {'rho_J':>8} {'linelet':>8} {'rho_L':>8}")
for ar in (1, 10, 100, 1000):
    nj, rj = count_sweeps(ar, 'jacobi')
    nl, rl = count_sweeps(ar, 'linelet')
    print(f"{ar:>6} {nj:>8} {rj:>8.4f} {nl:>8} {rl:>8.4f}")

Esto es lo que imprime.

    AR   Jacobi    rho_J  linelet    rho_L
     1       18   0.4870        8   0.1839
    10      514   0.9780        7   0.1699
   100     4879   0.9975        4   0.0194
  1000     5471   0.9978        2   0.0002

Al reducir $\Delta y_{\text{wall}}$ las cadenas se hunden más en la capa. Al subir la razón de crecimiento la malla vuelve antes a ser isótropa y las cadenas se acortan. Las cadenas son cortas. No cubren ni la mitad de las capas. El linelet sale barato por esa longitud.

Al aplastar mil veces, Jacobi se vuelve 300 veces más lento y choca contra el techo cercano a 6.000 que predijo la sección anterior. El linelet, en cambio, se acelera. El acoplamiento entre cadenas vecinas decae como $a_E/a_N = \mathrm{AR}^{-2}$, de modo que las cadenas se independizan entre sí. La cota del peor modo, $2a_E/(V/\Delta t + 2a_E) = 1/3$, sigue cumpliéndose; el decaimiento medido simplemente es mucho mejor.

Conviene partir solo la columna de Jacobi. 18, 514, 4879, 5471. Las relaciones de aspecto son 1, 10, 100 y 1000.

El primer factor diez en la relación multiplica el trabajo por 28,6. El siguiente factor diez lo multiplica por 9,5. El último lo multiplica por 1,12, que es 5471 dividido entre 4879. El cociente de coeficientes sigue creciendo como $\mathrm{AR}^2$ y el número de barridos deja de seguirlo. Con $N_y = 48$ la fórmula modal $1 - \frac{\pi^{2}}{2 N_y^{2}}$ vale 0,99786. Al resolver $0{,}99786^{n} = 10^{-6}$, n queda cerca de 6400. El 5471 impreso está justo debajo de ese techo. Aunque la malla se aplaste más, esta columna se queda pegada ahí.

La columna del linelet va al revés. 8, 7, 4, 2. Los factores medidos son 0,1839, 0,1699, 0,0194 y 0,0002. La cota $2a_E/(V/\Delta t + 2a_E) = 1/3$ no contiene $\mathrm{AR}$, y todo valor medido queda por debajo. En AR 1000, 0,0002 es cerca del 0,06% de 1/3. El acoplamiento entre cadenas, $a_E/a_N = \mathrm{AR}^{-2}$, se apaga, así que un paso de Thomas deja casi nada del error normal a la pared. En AR 1 son 8 barridos contra 18. La ganancia es pequeña. En AR 10 son 514 contra 7, unas 73 veces. Ahí se parte la elección.

La cota del suavizador puntual usa la misma diagonal. Con $V/\Delta t = 4a_E$ el denominador es $6a_E + 2a_N$. En AR 1 la cota vale 0,5 y el rho_J de la tabla es 0,4870. Casi el mismo número. En AR 100 la cota es 0,9998 y el rho_J medido es 0,9975. Cuando la cota se pega a 1, el 0,99786 de la fórmula modal es el número honesto. La estimación anterior de unos 6.000 barridos y el 5471 impreso son el mismo techo.

Trampas que conviene evitar al programarlo#

Tender las líneas al revés. Las cadenas deben seguir la dirección de coeficiente mayor, la normal a la pared. Tendidas a lo largo de la pared no aportan nada y solo cuestan el Thomas. Es fácil equivocarse, porque visualmente la malla parece "larga" en la dirección longitudinal.

Dejar $\sigma_{\min}$ en su valor por defecto. Si se cambia la malla sin revisar el umbral, las cadenas pueden desaparecer por completo. El precondicionador degenera entonces en Jacobi sin hacer ruido, y en el registro no aparece ninguna advertencia. Imprimir el número de linelets y su longitud media en cada ejecución detecta esto en una línea.

Dejar que las particiones MPI corten las cadenas. Una cadena seccionada en el borde de una partición queda en fragmentos, y la convergencia empeora cuantos más núcleos se añaden. Conviene dar al particionador un peso en la dirección de las líneas, o restringirlo para que no las corte.

Tapar un problema no lineal con el precondicionador. Si con linelet la simulación sigue divergiendo, el culpable suele ser el bucle externo y no el solver lineal. Primero hay que bajar el CFL y llevar el factor de subrelajación hacia 0,7 para reducir las propias actualizaciones; después se juzga. Un precondicionador resuelve más rápido la matriz que se le entrega, pero no arregla una matriz equivocada.

Thomas se desperdicia si las líneas no son normales a la pared#

Las cadenas siguen el coeficiente grande, solo en la normal a la pared. Tenderlas en la dirección del flujo es seguir con la vista el lado largo. Thomas sigue costando. El amortiguamiento coincide con Jacobi.

Conviene dejar en el registro el número de cadenas y la longitud media. Un cero impreso significa que $\sigma_{\min}$ ya no cuadra con la malla. El precondicionador ya es un suavizador diagonal. Sumar núcleos en ese estado solo parte más los fragmentos. El weak scaling que cae a un 16,3% de eficiencia mientras se refina la malla es la factura de ese corte.