---
title: Adjoints and sensitivities
summary: Gyre computes exact derivatives of outputs with respect to flow parameters and surface shape with one transposed linear solve per output. All partial derivatives come from automatic differentiation of the discrete operators, so gradients are consistent with the solver to solver tolerance.
section: theory
order: 8
level: 2
keywords: [discrete adjoint, total derivative, shape sensitivity, coupled adjoint, direct method, MDO]
---
<!-- gyre docs: expanded Markdown of theory/adjoint-and-sensitivities.md (generated tables, verification numbers and chart data are written out). Base level of this page: L2 Practice; passages at another level are bracketed by 'begin level N' / 'end level N' comments and start with [LN Name]. Built from revision ceb1f4b0be. -->

An optimizer needs $\mathrm dJ/\mathrm d\mathbf x$ for outputs $J$ (lift, drag, moment, loads) and
inputs $\mathbf x$ (angle of attack, Mach number, Reynolds number, the position of every surface
node). The state $\mathbf U$ is tied to $\mathbf x$ implicitly by $\mathbf R_w(\mathbf U, \mathbf x) = 0$.

<!-- begin level 2 (Practice) -->
**[L2 Practice] Two ways to get the same number**

$$
\frac{\mathrm dJ}{\mathrm d\mathbf x}
 = \frac{\partial J}{\partial \mathbf x} - \boldsymbol\psi^{\mathsf T}\frac{\partial \mathbf R_w}{\partial \mathbf x},
\qquad
\left(\frac{\partial \mathbf R_w}{\partial \mathbf U}\right)^{\!\mathsf T}\boldsymbol\psi = \left(\frac{\partial J}{\partial \mathbf U}\right)^{\!\mathsf T}
\qquad\text{(adjoint)}
$$

$$
\frac{\mathrm dJ}{\mathrm dx_k}
 = \frac{\partial J}{\partial x_k} + \frac{\partial J}{\partial \mathbf U}\frac{\mathrm d\mathbf U}{\mathrm dx_k},
\qquad
\frac{\partial \mathbf R_w}{\partial \mathbf U}\frac{\mathrm d\mathbf U}{\mathrm dx_k} = -\frac{\partial \mathbf R_w}{\partial x_k}
\qquad\text{(direct)}
$$

The adjoint costs one linear solve per *output* and is independent of the number of inputs —
the only practical choice for shape design, where inputs number in the hundreds to millions.
The direct method costs one solve per *input*. `totals(mode="auto")` picks the cheaper and both
are verified against each other and against finite differences.

**Table: Adjoint, finite-difference and direct derivatives** (Adjoint gradients versus finite differences)

| case | derivative | adjoint | finite difference | rel. diff. (adjoint, FD) | rel. diff. (adjoint, direct) |
|---|---|---|---|---|---|
| Euler NACA 0012 | dcl/dalpha | 7.59626 | 7.59626 | 9.8180e-11 | 2.7336e-09 |
| Euler NACA 0012 | dcd/dalpha | 0.0222355 | 0.0222355 | 2.0356e-08 | 2.8444e-08 |
| Euler NACA 0012 | dcm/dalpha | 0.0123702 | 0.0123702 | 1.1288e-08 | 4.6115e-07 |
| Euler NACA 0012 | dcl/dmach | 0.218262 | 0.218262 | 1.9822e-09 | 2.6393e-06 |
| Euler NACA 0012 | dcd/dmach | 3.8734e-04 | 3.8734e-04 | 2.1349e-08 | 1.1047e-05 |
| Euler NACA 0012 | dcm/dmach | 0.00128519 | 0.00128519 | 7.8285e-09 | 6.9061e-05 |
| Euler NACA 0012 | dcl/dshape (directional) | 0.0130491 | 0.0130491 | 6.7282e-09 | — |
| Euler NACA 0012 | dcd/dshape (directional) | -2.4501e-05 | -2.4501e-05 | 1.6480e-08 | — |
| Euler NACA 0012 | dcm/dshape (directional) | 3.9401e-04 | 3.9401e-04 | 3.8413e-08 | — |
| Cylinder Re = 40 | dcd/dre | -0.0125001 | -0.0125001 | 4.7607e-10 | — |
| Cylinder Re = 40 | dcd/dmach | 0.270411 | 0.270411 | 1.3254e-09 | — |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T11:18:25-0500; device cpu; quick suite. **Revision check: Same code as this build.** Produced by 39979fcf902ad6116eb550b881b4278d4599c8e2. Every commit between it and the commit this site was built from (ceb1f4b0bed41de53e586df10858079f5a2bcdd8) changes only documentation text, evidence records or bookkeeping: the code is the same. Full record: [adjoint_gradients.json](../data/verification/adjoint_gradients.json).

<!-- end level 2 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Where the partial derivatives come from**

Nothing is hand-differentiated. With $\mathbf d$ the surface displacement field and $\theta$ the
flow parameters, the code defines two pure functions,

```python
R(U, d, th) = residual_weak(U, geometry(deform(d)), th)
J(U, d, th) = functional(U, geometry(deform(d)), th)
```

and obtains

* $\partial J/\partial\mathbf U$ by `jax.grad`;
* $\boldsymbol\psi$ by FGMRES on the transposed operator from `jax.vjp(R, U)`, preconditioned by the
  transposed p-multigrid hierarchy;
* all remaining terms in a *single* reverse sweep, the gradient of the Lagrangian
  $\mathcal L(\mathbf d, \theta) = J(\mathbf U, \mathbf d, \theta) - \boldsymbol\psi^{\mathsf T}\mathbf R(\mathbf U, \mathbf d, \theta)$
  with $\mathbf U$ and $\boldsymbol\psi$ held fixed.

The last step delivers the derivative with respect to *every* surface node and every flow
parameter at the cost of roughly one residual evaluation. Shape sensitivities pass through mesh
deformation, metric terms, normals, boundary conditions and the functional's own geometric
dependence automatically.

<!-- end level 3 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Accuracy: exact for the discrete problem**

The adjoint derivative differs from the true derivative of the discrete output only through the
tolerances of the nonlinear and linear solves. Finite differences, by contrast, have a
truncation error that falls like $\varepsilon^2$ and a round-off error that grows like
$1/\varepsilon$; they agree with the adjoint only in a narrow window of step sizes:

**Chart: Finite-difference error of dcl/dα versus step size** (Adjoint gradients versus finite differences). Line chart of |FD − adjoint| / |adjoint| against step size, x axis log, y axis log; 1 series. The table holds the plotted data. Reference slopes drawn: 2 (ε² truncation). Smallest difference in this run: 1.2e-10 at a step of 1e-05. The reference line is the ε² truncation error of a central difference; round-off of the re-converged flows grows as 1/ε.

| step size | central difference |
|---|---|
| 0.01 | 1.0052e-04 |
| 0.001 | 9.7679e-07 |
| 1.0000e-04 | 9.6958e-09 |
| 1.0000e-05 | 1.1864e-10 |
| 1.0000e-06 | 1.8974e-10 |
| 1.0000e-07 | 2.3724e-10 |
| 1.0000e-08 | 2.7589e-09 |
| 1.0000e-09 | 3.9518e-07 |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T11:18:25-0500; device cpu; quick suite. **Revision check: Same code as this build.** Produced by 39979fcf902ad6116eb550b881b4278d4599c8e2. Every commit between it and the commit this site was built from (ceb1f4b0bed41de53e586df10858079f5a2bcdd8) changes only documentation text, evidence records or bookkeeping: the code is the same. Full record: [adjoint_gradients.json](../data/verification/adjoint_gradients.json).

This matters in optimization: gradient noise from badly chosen steps or loosely converged
solves is a classical cause of stalled optimizers.

<!-- end level 3 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Shape as an input: surface nodes, deformation, loads**

The shape input `x_surface` is the displacement of the high-order geometry nodes on the wall
(shape `(n_surface_nodes, dim)`); `surface_coordinates()` returns their baseline positions. A
geometry engine or parameterisation outside Gyre maps design variables to these displacements
and applies the chain rule to the returned sensitivity field.

Volume nodes follow by inverse-distance weighting over all boundary nodes: moving-wall nodes
carry their displacement, all other boundary nodes carry zero. The map is linear,
$\mathbf X = \mathbf X_0 + \mathbf W\mathbf d$, with a weight matrix that depends on the
baseline mesh only. $\mathbf W$ is dense, volume nodes by surface nodes: small meshes store
it, and above $2\times10^7$ entries it is applied matrix-free, recomputing the weights of a
chunk of nodes each time. The transposed map $\mathbf W\T$, which every shape derivative
needs, is written out the same way, because reverse-mode differentiation through the chunked
loop would keep the weights of every chunk, that is, the dense matrix again. Both
representations are the same map to round-off, in value, tangent and transpose.
On periodic meshes the weighting is periodic: distances are taken to the nearest periodic
image, twin nodes on the two periodic planes are one source and share one row of $\mathbf W$,
so the two sides of a periodic interface stay matched for any surface displacement. Twin
surface nodes keep their own rows of $\mathbf d$ but move by the mean of the two prescribed
values, and the derivative with respect to $\mathbf d$ is shared equally between them. Only
translational periodicity is supported.

The matching output for structural coupling is `f_surface`: consistent nodal loads obtained by
integrating the wall traction against the geometry shape functions. They sum exactly to the
total force and are work-conjugate to `x_surface`.

<!-- end level 3 -->

<!-- begin level 4 (Theory) -->
**[L4 Theory] Building a coupled adjoint from reduced products**

For a coupled aero-structural system $\mathbf d = \mathbf C\,\mathbf f(\mathbf d, \alpha)$ with objective
$J(\mathbf d, \alpha)$ the total derivative is

$$
\frac{\mathrm dJ}{\mathrm d\alpha}
 = \frac{\partial J}{\partial\alpha} + \boldsymbol\lambda^{\mathsf T}\mathbf C\,\frac{\partial\mathbf f}{\partial\alpha},
\qquad
\boldsymbol\lambda = \frac{\partial J}{\partial\mathbf d}^{\mathsf T} + \left(\frac{\partial\mathbf f}{\partial\mathbf d}\right)^{\!\mathsf T}\mathbf C^{\mathsf T}\boldsymbol\lambda ,
$$

where all flow derivatives are *total* derivatives at converged flow (the flow state is
eliminated). The flow discipline has to supply exactly two operations:

* `totals(of=J, wrt=[alpha, x_surface])` — the reduced gradient of a scalar;
* `adjoint_product({"f_surface": seed}, wrt=[...])` — the reverse product
  $(\mathrm d\mathbf f/\mathrm d\mathbf x)^{\mathsf T}\mathbf{seed}$, one flow-adjoint solve whatever the
  number of load components.

The second is what `FlowDiscipline.adjoint_product` provides (and its forward counterpart
`tangent_product`). The [aero-structural example](../examples/aerostructural-coupling.md)
iterates the $\boldsymbol\lambda$ equation with these calls and matches finite differences of the
fully re-converged coupled system to about $10^{-6}$. Frameworks that prefer to own the linear
solve instead use the *partial* products `jvp`/`vjp` and `solve_linear`.

<!-- end level 4 -->

<!-- begin level 4 (Theory) -->
**[L4 Theory] Why the discrete adjoint, and why adjoint consistency still matters**

Differentiating the discrete equations (rather than discretizing the continuous adjoint
equations) guarantees that the gradient is the gradient of the function the optimizer actually
evaluates, at any resolution. It does not by itself guarantee that $\boldsymbol\psi$ approximates the
continuous adjoint: that additionally requires an [adjoint-consistent](boundary-conditions.md)
discretization. Gyre has both properties, which is why the same $\boldsymbol\psi$ can be reused as a
weight in [output error estimation](error-estimation.md).

<!-- end level 4 -->