[Paper Review] Tired of Differentiating Jacobians by Hand? — Operator-Overloading Automatic Differentiation (ADOO)

I started the HLLC flux Jacobian by hand and threw out three sheets. One product-rule term dropped, and the code diverged in silence. The flux was fine. The derivative was not. Half of implicit CFD is filling that matrix — the Jacobian of the spatial discretization — without a sign error. ADOO, operator-overloading automatic differentiation from Fraysse et al. (2019), hands the derivative to the code. A scheme with a root finder inside, Godunov included, runs without a hand formula. How much Newton gains from that is written on the convergence curve.

Shrink the finite-difference step h and the error curve dips once, then climbs again. Drag x and the dual path stays pinned to the floor.

The bottom of the valley sits near $h\sim 10^{-8}$. Leave that pocket and only the amber curve climbs back up.

Paper Information#

• Title: Automatic Differentiation using Operator Overloading (ADOO) for implicit resolution of hyperbolic single phase and two-phase flow models
• Authors: G. Fraysse et al.
• Year: 2019
• Keywords: automatic differentiation, implicit, two-phase, finite volume, unstructured meshes

One-line summary: leave the flux code untouched, swap the data type, and read off the exact Jacobian.

Why the Jacobian Is a Headache#

Implicit time stepping solves a nonlinear residual $\mathbf{R}(\mathbf{Q}) = 0$ at each step with Newton iterations.

Here $\mathbf{J}$ is the residual Jacobian (the matrix of partials with respect to each conserved variable) and $\delta\mathbf{Q}$ is the update. For Newton to converge quadratically, $\mathbf{J}$ has to be exact.

There are three routes, and each has a catch. Hand differentiation (analytic) is exact but has to be re-derived whenever you touch the scheme. For a branch-heavy flux like AUSM+ or an iterative one like Godunov, the derivation is a nightmare. Finite differences reuse the code but leave you stuck choosing a step $h$ between round-off and truncation error.

Newton-Krylov-matrix-free only finite-differences the Jacobian-vector product, so it never builds the matrix. Stopping a Krylov run that otherwise keeps going for a thousand steps still wants the actual matrix entries. ADOO cuts the knot by computing the derivative instead of approximating it.

Dual Numbers: Carrying a Derivative Alongside the Value#

The idea is plain. Replace a real number $x$ with an object of two components $(v, dv)$: $v$ is the value, $dv$ is its derivative at that point. That is a dual number.

The arithmetic rules are just the product and chain rules, transcribed.

The left of each rule is the value; the right updates the derivative component. Take $f(x) = \sin(x^2) + x^2$. By hand, $f'(x) = \cos(x^2)\cdot 2x + 2x$. With a dual number, set $x$'s derivative component to $1$ (since $dx/dx = 1$), run the code, and the final object's $dv$ is $f'(x)$.

from dataclasses import dataclass
import math
 
@dataclass
class Dual:
    v: float   # value
    d: float   # derivative component
 
    # operator overloading — each rule is product / chain rule (paper 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)
 
# differentiate f(x) = sin(x^2) + x^2 at x = 1.3
x = Dual(1.3, 1.0)          # seed: dx/dx = 1
f = sin_d(x * x) + x * x    # the original expression, unchanged
print(f.v, f.d)             # value, f'(1.3)
print(math.cos(1.3**2) * 2 * 1.3 + 2 * 1.3)  # cross-check vs hand derivative

f.d matches the hand derivative down to round-off. The hand-derived formula here is only a check — it never entered the actual computation. That's the whole point.

Why You Can't Trust Finite Differences#

The real value of ADOO shows up next to finite differences. The error of a central difference $[f(x+h)-f(x-h)]/2h$ is a tug of war. Large $h$ lets truncation error (the Taylor remainder) dominate; too-small $h$ makes round-off (loss of significant digits when subtracting two nearby numbers) explode. So the error traces a U shape in $h$. The sweet spot sits near $h \sim \sqrt{\epsilon} \approx 10^{-8}$, but even that shifts from problem to problem.

The amber curve above is that U, and the cyan dashed line is the floor with no step on it. AD has no step, so there is nothing to tune.

The Euler Flux Jacobian: By Hand vs AD#

Now leave the scalar behind. The 1D compressible Euler system has conserved variables $\mathbf{Q} = [\rho, \rho u, \rho E]^\top$ and flux $\mathbf{F}(\mathbf{Q})$. For an ideal gas the flux Jacobian $\partial\mathbf{F}/\partial\mathbf{Q}$ is the well-known $3\times3$ matrix, but hand-deriving it tangles $\gamma$ with kinetic-energy terms. Extend the dual number to a vector (make each $dv$ a 3-component array) and one run of the code fills the whole matrix.

import numpy as np
 
class DualVec:
    """value + gradient (seed vector) w.r.t. 3 independent variables."""
    def __init__(self, v, grad):
        self.v = v
        self.g = np.asarray(grad, float)  # d(this)/dQ, length 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)  # product rule
        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):
    # seed each conserved variable: 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                                  # velocity
    kinetic = rhou * u * 0.5                         # ½ρu²
    p = (rhoE - kinetic) * (gamma - 1.0)             # pressure (ideal gas)
 
    F0 = rhou                                        # ρu
    F1 = rhou * u + p                                # ρu² + p
    F2 = (rhoE + p) * u                              # (ρE + p)u
 
    # the .g of each flux component IS that row of the Jacobian (paper eqs 9–10)
    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)
 
# cross-check against the hand-derived exact Jacobian
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("max error:", np.abs(J_ad - J_ref).max())   # ~1e-15

euler_flux_jacobian follows the same arithmetic the Riemann solver's flux code actually uses. We never derived the Jacobian, yet its accuracy is machine precision. Swap the scheme for AUSM+? Replace only the flux function and the Jacobian follows for free.

Plug the code's state into the hand formula and the matrix closes in numbers. With $Q=(1.2,\ 0.6,\ 3.0)$ and $\gamma=1.4$, $u=0.6/1.2=0.5$. The momentum row is $0.5(\gamma-3)u^2$, $(3-\gamma)u$, $\gamma-1$, which is $(-0.2,\ 0.8,\ 0.4)$. The energy row is longer. $(\gamma-1)u^3=0.05$ and $\gamma u E/\rho=1.75$, so the first entry is $0.05-1.75=-1.70$. From $\gamma E/\rho=3.5$, subtract $1.5(\gamma-1)u^2=0.15$ and the middle entry is $3.35$. The last entry is $\gamma u=0.7$. The mass row stays $(0,\ 1,\ 0)$. The code comment puts the max gap between this matrix and AD at $\sim 10^{-15}$. The scalar check sits in the same place. At $x=1.3$, $x^2=1.69$ and $2x=2.6$, so the hand derivative is $\cos(1.69)\cdot 2.6+2.6$, and f.d splits from that only inside round-off.

Newton Convergence: What Exactness Buys#

Why the exact Jacobian matters is written in the convergence curve. An approximate Jacobian drops Newton from quadratic to linear convergence. Iteration counts climb, and at large CFL it just diverges. The paper reports that with an exact Jacobian it reaches a $10^{-6}$ residual in fewer than 10 iterations even at CFL 20 — the signature of quadratic convergence.

Below, push the Jacobian-error slider up from zero.

At 0% error (the AD-exact Jacobian) the residual roughly doubles its digits per iteration and plunges — the signature of quadratic convergence. Add just 20% error and the curve flattens into a gentle straight line (first order), taking many more iterations to reach the same precision. Factor in the Krylov iterations and that gap becomes wall-clock time.

A Critical Look: The Shadow Side of ADOO#

ADOO isn't free. Operator overloading builds an object per scalar and multiplies arrays. Fraysse's paper admits that forward-mode cost scales with the number of independent variables — for a large-block multiphase system those arrays get heavy. That's exactly why source-code transformation (ADSCT, e.g. Tapenade) can be faster through compile-time optimization.

I hit a practical issue while reproducing this. When you differentiate a scheme with a root-finding loop, like Godunov, you have to converge not just the value but the derivative component too. The derivative usually converges slower than the value, so a stopping criterion based on the value alone leaves the Jacobian contaminated. That one sentence in the paper saved me days.

From an OpenFOAM standpoint this isn't alien. There have been experiments assembling blockLduMatrix with a dual scalar type, and SU2 already builds its adjoint with code-transformation AD. The message is the same — the era of humans deriving Jacobians is fading.

Stopping the root finder on the value alone dirties the Jacobian#

A root finder inside the flux, Godunov being the usual case, does not converge the value and the derivative at the same rate. The derivative is slower. Stop on the value and the flux looks right while the matrix entries are already dirty. The paper's exact Jacobian still reaches a residual of $10^{-6}$ in under 10 iterations at CFL 20. Mix in a 20% error and quadratic convergence lies down into a straight line. Digits that were roughly doubling each iteration become a constant slope.

Forward-mode cost scales with the number of independent variables. Large multiphase blocks make the dual arrays heavy, and source transformation can be the cheaper route. A branched flux, the kind used to fill an interface with a ghost fluid, is where hand differentiation dies first. Leave the stopping test on the value alone and that Jacobian walks into Newton already wrong.