[Reseña] ¿Cansado de derivar jacobianos a mano? — Diferenciación automática por sobrecarga de operadores (ADOO)

La derivación a mano del jacobiano del flujo HLLC acabó en tres hojas tiradas. Cayó un término de la regla del producto y el código divergió en silencio. El flujo estaba bien. La derivada no. La mitad del CFD implícito consiste en llenar esa matriz —el jacobiano de la discretización espacial— sin un signo erróneo. ADOO, la diferenciación automática por sobrecarga de operadores de Fraysse et al. (2019), deja la derivada en el código. Un esquema con una búsqueda de raíces adentro, Godunov incluido, corre sin fórmula manual. Cuánto gana Newton con eso lo escribe la curva de convergencia.

Al reducir el paso h de las diferencias finitas, la curva de error baja una vez y vuelve a subir. Al mover x, la ruta dual sigue pegada al piso.

El fondo del valle queda cerca de $h\sim 10^{-8}$. Fuera de ese bolsillo, solo la curva ámbar vuelve a trepar.

Información del artículo#

• Título: Automatic Differentiation using Operator Overloading (ADOO) for implicit resolution of hyperbolic single phase and two-phase flow models
• Autores: G. Fraysse et al.
• Año: 2019
• Palabras clave: automatic differentiation, implicit, two-phase, finite volume, unstructured meshes

Resumen en una línea: dejar intacto el código del flujo, cambiar solo el tipo de dato y leer el jacobiano exacto.

Por qué el jacobiano es un dolor de cabeza#

El avance temporal implícito resuelve un residuo no lineal $\mathbf{R}(\mathbf{Q}) = 0$ en cada paso con iteraciones de Newton.

Aquí $\mathbf{J}$ es el jacobiano del residuo (la matriz de derivadas parciales respecto a cada variable conservada) y $\delta\mathbf{Q}$ es la actualización. Para que Newton converja de forma cuadrática, $\mathbf{J}$ tiene que ser exacto.

Hay tres caminos, y cada uno tiene su trampa. La derivación a mano (analítica) es exacta, pero hay que rederivarla cada vez que se toca el esquema. Para un flujo lleno de bifurcaciones como AUSM+ o uno iterativo como Godunov, la derivación misma es una pesadilla. Las diferencias finitas reutilizan el código, pero obligan a elegir un paso $h$ entre el error de redondeo y el de truncamiento.

Newton-Krylov-matrix-free solo aproxima por diferencias finitas el producto jacobiano-vector, así que nunca arma la matriz. Para detener un Krylov que, sin matriz de precondicionamiento, sigue mil pasos siguen haciendo falta las entradas reales. ADOO corta el nudo calculando la derivada en lugar de aproximarla.

Números duales: llevar la derivada junto al valor#

La idea es sencilla. Se reemplaza un número real $x$ por un objeto de dos componentes $(v, dv)$: $v$ es el valor, $dv$ es su derivada en ese punto. Eso es un número dual.

Las reglas aritméticas son simplemente la regla del producto y la de la cadena, transcritas.

El lado izquierdo de cada regla es el valor; el derecho actualiza la componente de la derivada. Tomemos $f(x) = \sin(x^2) + x^2$. A mano, $f'(x) = \cos(x^2)\cdot 2x + 2x$. Con un número dual, se fija la componente de derivada de $x$ en $1$ (ya que $dx/dx = 1$), se ejecuta el código, y el $dv$ del objeto final es $f'(x)$.

from dataclasses import dataclass
import math
 
@dataclass
class Dual:
    v: float   # valor
    d: float   # componente de derivada
 
    # sobrecarga de operadores — cada regla es la del producto / cadena (artículo III.2)
    def __add__(self, o):
        o = o if isinstance(o, Dual) else Dual(o, 0.0)
        return Dual(self.v + o.v, self.d + o.d)
 
    def __mul__(self, o):
        o = o if isinstance(o, Dual) else Dual(o, 0.0)
        return Dual(self.v * o.v, self.v * o.d + self.d * o.v)
 
    def __truediv__(self, o):
        o = o if isinstance(o, Dual) else Dual(o, 0.0)
        return Dual(self.v / o.v, (self.d * o.v - self.v * o.d) / (o.v * o.v))
 
def sin_d(a: Dual) -> Dual:
    return Dual(math.sin(a.v), math.cos(a.v) * a.d)
 
# derivar f(x) = sin(x^2) + x^2 en x = 1.3
x = Dual(1.3, 1.0)          # semilla: dx/dx = 1
f = sin_d(x * x) + x * x    # la expresión original, sin tocar
print(f.v, f.d)             # valor, f'(1.3)
print(math.cos(1.3**2) * 2 * 1.3 + 2 * 1.3)  # contraste con la derivada a mano

El f.d de la salida coincide con la derivada a mano hasta el nivel del redondeo. La fórmula derivada a mano aquí es solo una comprobación: nunca entró en el cálculo real. Ese es todo el punto.

Por qué no se puede confiar en las diferencias finitas#

El verdadero valor de ADOO aparece junto a las diferencias finitas. El error de una diferencia central $[f(x+h)-f(x-h)]/2h$ es un tira y afloja. Un $h$ grande deja que domine el error de truncamiento (el resto de Taylor); un $h$ demasiado pequeño hace estallar el redondeo (pérdida de dígitos significativos al restar dos números cercanos). Así, el error dibuja una U en $h$. El punto óptimo está cerca de $h \sim \sqrt{\epsilon} \approx 10^{-8}$, pero hasta eso se corre según el problema.

La curva ámbar de arriba es esa U, y la línea punteada cian es el piso sin paso. AD no tiene paso, así que no hay nada que ajustar.

El jacobiano del flujo de Euler: a mano vs AD#

Ahora dejemos el escalar. El Euler compresible 1D tiene variables conservadas $\mathbf{Q} = [\rho, \rho u, \rho E]^\top$ y flujo $\mathbf{F}(\mathbf{Q})$. Para un gas ideal, el jacobiano del flujo $\partial\mathbf{F}/\partial\mathbf{Q}$ es la conocida matriz $3\times3$, pero derivarla a mano enreda $\gamma$ con los términos de energía cinética. Al extender el número dual a un vector (haciendo que cada $dv$ sea un arreglo de 3 componentes), una sola ejecución del código llena la matriz entera.

import numpy as np
 
class DualVec:
    """valor + gradiente (vector semilla) respecto a 3 variables independientes."""
    def __init__(self, v, grad):
        self.v = v
        self.g = np.asarray(grad, float)  # ∂(this)/∂Q, longitud 3
    def __add__(s, o):  return DualVec(s.v + o.v, s.g + o.g)
    def __sub__(s, o):  return DualVec(s.v - o.v, s.g - o.g)
    def __mul__(s, o):
        if isinstance(o, DualVec):
            return DualVec(s.v * o.v, s.v * o.g + s.g * o.v)  # regla del producto
        return DualVec(s.v * o, s.g * o)
    def __truediv__(s, o):
        return DualVec(s.v / o.v, (s.g * o.v - s.v * o.g) / (o.v * o.v))
 
def euler_flux_jacobian(Q, gamma=1.4):
    # sembrar cada variable conservada: rho -> (Q0, e0), etc.
    rho  = DualVec(Q[0], [1, 0, 0])
    rhou = DualVec(Q[1], [0, 1, 0])
    rhoE = DualVec(Q[2], [0, 0, 1])
 
    u = rhou / rho                                  # velocidad
    kinetic = rhou * u * 0.5                         # ½ρu²
    p = (rhoE - kinetic) * (gamma - 1.0)             # presión (gas ideal)
 
    F0 = rhou                                        # ρu
    F1 = rhou * u + p                                # ρu² + p
    F2 = (rhoE + p) * u                              # (ρE + p)u
 
    # el .g de cada componente de F ES esa fila del jacobiano (ecs. 9–10 del artículo)
    return np.array([F0.g, F1.g, F2.g])
 
Q = np.array([1.2, 0.6, 3.0])   # ρ, ρu, ρE
J_ad = euler_flux_jacobian(Q)
 
# contraste con el jacobiano exacto derivado a mano
rho, mom, E = Q
u = mom / rho; g = 1.4
J_ref = np.array([
    [0, 1, 0],
    [0.5*(g-3)*u*u, (3-g)*u, g-1],
    [((g-1)*u**3 - g*u*E/rho), (g*E/rho - 1.5*(g-1)*u*u), g*u],
])
print("error máximo:", np.abs(J_ad - J_ref).max())   # ~1e-15

euler_flux_jacobian sigue la misma aritmética que el código de flujo que el solucionador de Riemann usa de verdad. Nunca se derivó el jacobiano, y aun así su exactitud es de precisión de máquina. ¿Cambiar el esquema por AUSM+? Basta reemplazar la función de flujo y el jacobiano viene solo.

Al meter el estado del código en la fórmula manual, la matriz se cierra en números. Con $Q=(1.2,\ 0.6,\ 3.0)$ y $\gamma=1.4$, $u=0.6/1.2=0.5$. La fila de cantidad de movimiento es $0.5(\gamma-3)u^2$, $(3-\gamma)u$, $\gamma-1$, es decir $(-0.2,\ 0.8,\ 0.4)$. La fila de energía es más larga. $(\gamma-1)u^3=0.05$ y $\gamma u E/\rho=1.75$, así que la primera entrada vale $0.05-1.75=-1.70$. Desde $\gamma E/\rho=3.5$, al restar $1.5(\gamma-1)u^2=0.15$ la entrada central queda en $3.35$. La última es $\gamma u=0.7$. La fila de masa sigue en $(0,\ 1,\ 0)$. El comentario del código sitúa la diferencia máxima entre esta matriz y AD en $\sim 10^{-15}$. La comprobación escalar ocupa el mismo sitio. En $x=1.3$, $x^2=1.69$ y $2x=2.6$, de modo que la derivada manual es $\cos(1.69)\cdot 2.6+2.6$, y f.d solo se separa de ese valor dentro del redondeo.

Convergencia de Newton: lo que compra la exactitud#

Por qué importa el jacobiano exacto lo dice la curva de convergencia. Un jacobiano aproximado baja Newton de convergencia cuadrática a lineal. El número de iteraciones sube, y con CFL grande simplemente diverge. El artículo reporta que, gracias al jacobiano exacto, alcanza un residuo de $10^{-6}$ en menos de 10 iteraciones incluso a CFL 20: la firma de la convergencia cuadrática.

Al subir desde cero el control de error del jacobiano, cambia el orden de la curva.

Con 0% de error (el jacobiano exacto de AD), el residuo duplica más o menos sus dígitos por iteración y se desploma: la firma de la convergencia cuadrática. Con solo 20% de error, la curva se acuesta en una recta suave (primer orden) y hacen falta muchas más iteraciones para la misma precisión. Si se cuentan las iteraciones de Krylov, esa brecha se vuelve tiempo de reloj.

Una mirada crítica: la sombra de ADOO#

ADOO no es gratis. La sobrecarga de operadores crea un objeto por escalar y multiplica arreglos. El artículo de Fraysse admite que el costo del modo directo (forward) crece con el número de variables independientes: en un sistema multifásico de bloques grandes esos arreglos pesan. Por eso la transformación de código fuente (ADSCT, p. ej. Tapenade) puede ser más rápida gracias a la optimización en tiempo de compilación.

Al reproducirlo surgió un problema práctico. Al diferenciar un esquema con un bucle de búsqueda de raíces, como Godunov, hay que converger no solo el valor sino también la componente de la derivada. La derivada suele converger más lento que el valor, así que un criterio de parada basado solo en el valor deja el jacobiano contaminado. Esa frase del artículo me ahorró días.

Desde la óptica de OpenFOAM esto no es ajeno. Hubo experimentos ensamblando blockLduMatrix con un tipo escalar dual, y SU2 ya construye su adjunto con AD por transformación de código. El mensaje es el mismo: la era de los jacobianos derivados por humanos se apaga.

Parar la búsqueda de raíces solo con el valor ensucia el jacobiano#

Una búsqueda de raíces dentro del flujo, el caso habitual de Godunov, no converge el valor y la derivada al mismo ritmo. La derivada va más lenta. Si el corte mira solo el valor, el flujo parece correcto y las entradas de la matriz ya están sucias. El jacobiano exacto del artículo todavía llega a un residuo de $10^{-6}$ en menos de 10 iteraciones con CFL 20. Mezclar un 20% de error tumba la convergencia cuadrática en una recta. Los dígitos que casi se duplicaban por iteración pasan a una pendiente casi constante.

El costo del modo directo crece con el número de variables independientes. Los bloques multifásicos grandes vuelven pesados los arreglos duales, y la transformación de código puede salir más barata. Un flujo con muchas ramas, del tipo que rellena una interfaz con un fluido fantasma, es donde la derivación manual muere primero. Si el criterio de paro no llega hasta la componente de la derivada, ese jacobiano entra en Newton ya equivocado.