Integrating with an MDO framework

1. In-process (Python). Construct a case and use its discipline directly:

from gyre.cases import Case, preset

flow = Case(preset("naca0012_euler")).discipline
flow.set_inputs(alpha=0.04, mach=0.55)
flow.solve()
out = flow.outputs(("cl", "cd", "f_surface"))
grad = flow.totals(of=("cl", "cd"), wrt=("alpha", "x_surface"))["derivatives"]

2. Remote session (HTTP). The framework runs on another machine; a warm solver process lives in a Gyre deployment:

from gyre.client import GyreClient

with GyreClient("http://gyre-host:8600").session("naca0012_euler") as flow:
    out = flow.evaluate(alpha_deg=2.5)["outputs"]
    grad = flow.gradients(of=["cl"], wrt=["alpha_deg", "x_surface"])
    back = flow.adjoint({"f_surface": seed}, wrt=["x_surface", "alpha"])

The calls a driver makes on such a session — at fixed design points, without an optimizer — are shown by examples/remote_session.py, described under A remote session called at fixed design points.

RemoteDiscipline exposes the reduced-space contract under its own names: evaluate(**inputs) (set inputs, solve, outputs), gradients(of, wrt) (the derivatives of totals), adjoint(seeds, wrt) and tangent(d_inputs, outputs) (the reverse and forward total products), residual(**inputs) (norms of the residual of the session’s state), revert() after a failed evaluation, get_state / set_state, checkpoint() / restore() (a checkpoint also rebuilds a session that was closed or reaped: GyreClient.restore_session), plus describe, surface_coordinates, error_estimate and distributions. The partial products jvp, vjp and solve_linear are available in-process only. Fields travel as JSON arrays.

3. Framework adapter. gyre.integrations.openmdao wraps a discipline for OpenMDAO in two forms (see the OpenMDAO example): GyreFlowComp, a matrix-free ExplicitComponent (reduced space: the framework sees inputs and outputs), and GyreFlowGroup, an ImplicitComponent whose output is the flow state plus a functionals component (residual level: the framework sees the state, the residual and its linear solves). Both turn a flow solve that fails and a linear solve that stalls into an OpenMDAO analysis error.

• Non-dimensional free stream: $\rho_\infty = |\mathbf u_\infty| = T_\infty = 1$; coefficients are $2F/S_\text{ref}$ with $S_\text{ref} = 1$ in 2-D. Spanwise-extruded 3-D cases use $S_\text{ref} = $ span, so their coefficients, and the error estimates of those coefficients, are per unit span. reference.area, reference.length and reference.moment_center of the case set other reference quantities, which a case on a mesh file usually needs.
• alpha is in radians. Through the service alpha_deg is accepted as a session input (hence by RemoteDiscipline.evaluate) and in the wrt of session gradients and adjoint and of job gradients; the derivative is then per degree. The in-process FlowDiscipline knows only alpha.
• x_surface is a displacement relative to the baseline surface nodes returned by surface_coordinates(); f_surface is a force per node (not a traction), already integrated, and sums to the total force: 2 * f_surface.sum(0) == (cx, cy).
• Loads of components. reference.components of the case names sets of wall groups, for example {wing: [wing], tail: [htail, vtail]} on a mesh with those groups, and the discipline then returns every coefficient per component as well (cd_wing, cmy_tail, …) from the same solve; each can be an objective or a constraint with its own adjoint. reference.groups instead changes what the plain coefficients cover (default: all walls). f_surface and x_surface always cover every wall, so the identity above holds for the coefficients of all walls only.
• Derivatives with respect to x_surface have the same shape as x_surface.
• On a periodic mesh (spanwise-extruded cases) the surface nodes on the two periodic planes are one topological node. Both rows remain inputs, but the pair moves by the mean of the two prescribed displacements, and the derivative with respect to x_surface is shared equally between them. Give both rows the same value.

Two styles of multidisciplinary derivative computation exist, and Gyre supports both.

Reduced space (the flow state is eliminated inside the discipline). The framework sees an explicit map inputs → outputs with total derivative products. Each reverse product costs one flow-adjoint solve; this is what adjoint_product, tangent_product, the session endpoint /adjoint and the OpenMDAO adapter provide. It is the simplest correct integration and works across a network boundary because the state never leaves the solver.

Full space (the framework owns the coupled linear solve). The framework treats the flow as an implicit component and needs the residual, partial products of the residual and outputs, and a preconditioned linear solve, all about a state it owns: apply_nonlinear, apply_linear (jvp, vjp) and solve_linear with state= and inputs=. This is the MPhys-style contract; it is available in-process, and GyreFlowGroup is its OpenMDAO form.

For a structure $\mathbf d = \mathbf C\,\mathbf f$ coupled to the flow and an objective $J$:

tot = flow.totals(of=("cl",), wrt=("alpha", "x_surface"))["derivatives"]["cl"]
lam = tot["x_surface"].ravel()
for _ in range(max_iter):                       # transposed block Gauss–Seidel
    seed = (C.T @ lam).reshape(shape)
    back = flow.adjoint_product({"f_surface": seed}, wrt=("x_surface",))["derivatives"]
    lam = tot["x_surface"].ravel() + back["x_surface"].ravel()
cross = flow.adjoint_product({"f_surface": (C.T @ lam).reshape(shape)}, wrt=("alpha",))
dJ_dalpha = tot["alpha"] + cross["derivatives"]["alpha"]

The aero-structural example runs exactly this and verifies it against finite differences of the re-converged coupled system. The derivation is in Adjoints and sensitivities.

Gyre differentiates with respect to surface-node positions because that is the one interface every geometry engine can target. If your design variables $\boldsymbol\theta$ move the surface by $\mathbf d(\boldsymbol\theta)$ — free-form deformation, CAD parameters, Hicks–Henne bumps — then

one matrix–vector product with the field Gyre returns. Volume-mesh deformation is handled inside Gyre and is already part of $\mathrm dJ/\mathrm d\mathbf d$.

• Use sessions. The first evaluation of a mesh/order compiles kernels (tens of seconds); later evaluations reuse them and warm-start from the previous flow.
• Converge tightly when you need gradients. The adjoint is exact for the converged discrete problem; a loosely converged flow gives a gradient of something slightly different. Presets converge ten orders of magnitude by default.
• One adjoint per output. Requesting of=["cl","cd","cm"] costs three transposed solves; request only what the optimizer uses, or seed a linear combination with adjoint_product.
• On wall-resolved meshes, precondition the linear solves along the wall-normal lines. Every output of every design point costs one transposed solve to a tight tolerance, and solver.smoother: line is made for such solves (how, what it gained on the flat plate). The flow solve itself need not get faster with it, and no study measures how a design loop divides its time between flow solves and transposed solves: that depends on the number of outputs and on the step between design points. The smoother belongs to the case, so a session keeps it for every evaluation.
• Check the mesh for the objective, not in general. error_estimate tells you whether a change in the objective between two designs exceeds its discretization error.
• Failures are data. A non-converged evaluation returns converged: false with the reason and the outputs and residual of the iterate it ended at. The failed iterate is not accepted: the state stays the last converged solution and the next evaluation warm-starts from it. In Python, flow.U is unchanged, flow.history describes the failed attempt and flow.failed_state holds its last iterate. Repeating the same evaluation repeats the same warm start, so take a smaller step, or in Python start elsewhere with flow.solve(U0=...).
• No derivatives of a failed point. After a solve that did not converge the kept state belongs to the previous inputs, so totals, tangent_product, adjoint_product, jvp, vjp, solve_linear and linearize raise RuntimeError until a solve converges (or the caller declares a state with set_state). A session answers gradients, adjoint, error_estimate, distributions and the field request with 422 RuntimeError in the same situation, and GET /sessions/{id}/state reports converged: false. outputs() raises as well, because the kept state belongs to the previous inputs: name a state to get numbers anyway (outputs(state=flow.failed_state) is the failed iterate), or go back to the last converged point with revert_inputs() (revert() on a session).
• Choose the degree once, then keep it. adapt on a job, --target-error on the command line or gyre.adapt.adaptive_order select the polynomial degree from an error tolerance on one output (Gradients and error estimates). A session keeps the degree of its case.

All partial derivatives are obtained by automatic differentiation of the discrete operators, so there is no approximation in the Jacobian, the boundary conditions, the metric terms or the mesh deformation. Total derivatives are exact up to the tolerance of the linear solve (default relative $10^{-10}$) and the convergence of the nonlinear solve. They are derivatives of the discrete outputs — the function the optimizer actually evaluates — and therefore consistent with line searches and convergence checks at any resolution.