Adjoints and sensitivities

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.

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

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.

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:

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

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.

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

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 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.

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 discretization. Gyre has both properties, which is why the same $\boldsymbol\psi$ can be reused as a weight in output error estimation.