Design decisions

Context. An MDO discipline needs Jacobian-vector products for Newton–Krylov, transposed products for adjoints, and derivatives with respect to geometry and flow parameters. Hand-differentiated solvers carry a permanent maintenance burden and a permanent risk that the derivative code lags the primal code.

Decision. Write every numerical kernel as a pure function of JAX arrays; obtain all derivatives with jax.jvp, jax.vjp, jax.linearize and jax.grad; vectorise with array programming and vmap.

Cost. Just-in-time compilation: tens of seconds for each new combination of mesh size and polynomial degree, which dominates small runs. Mitigations: a persistent compilation cache, warm sessions, CPU execution for small problems. Static shapes rule out some adaptive algorithms (e.g. per-element polynomial degree) without restructuring.

Rejected. A C++ core with source-transformation or operator-overloading AD: faster cold start, but far more code, a separate adjoint build, and no free GPU path.

Context. High-order methods are attractive for design (accuracy per degree of freedom, smooth dependence of outputs on inputs) but classical collocation schemes are fragile when under-resolved.

Decision. Flux-differencing DGSEM on Gauss–Lobatto nodes with entropy-conservative two-point fluxes in the volume and upwind dissipation at interfaces.

Cost. Two-point fluxes are evaluated for node pairs rather than nodes (mitigated by the symmetric pair tables), and Gauss–Lobatto collocation is somewhat less accurate per degree of freedom than Gauss collocation.

Rejected. Standard (over-integrated) DG: robust only with expensive over-integration and no discrete entropy statement. Finite volume: far simpler, but low order.

Context. Matrix-free Krylov solvers need a preconditioner; the workhorse for DG is element block Jacobi, usually inside multigrid. Extracting blocks from an AD Jacobian normally needs graph colouring and many JVPs, or sparse assembly.

Decision. Give the residual two state arguments: one used for everything element-local, one for everything gathered from neighbours. The partial derivative with respect to the first is exactly block diagonal, so block_size JVPs with identical per-element unit tangents return all blocks.

Cost. Every boundary condition and every face term must respect the convention (covered by a test that compares against a dense Jacobian for each wall type), and each interior face is evaluated from both sides.

Context. A capability exposed through four interfaces tends to grow four subtly different behaviours; and a numerical solver inside a web server can stall or crash it.

Decision. All interfaces construct a CaseSpec and call the same runner. The API process never executes solver kernels: jobs run in fresh worker processes, sessions in long-lived ones, each with its own JAX runtime and device selection. MCP tools call the REST API through an in-process ASGI transport, so agents exercise the same code path as programs.

Cost. Process start-up and JIT per job (hence sessions and the result cache), and data crossing process boundaries as JSON.

Context. A RANS model brings stiff sources, a wall distance and a variable six orders smaller than the mean flow. Solvers often treat it apart: loosely coupled, frozen in the adjoint, with a wall distance that is not differentiated.

Decision. The model variable is appended to the state and its equation to the residual (formulation). Sources are integrated at Gauss points with the BR2 gradient; the wall distance is stored as foot-point weights on mesh nodes and re-evaluated inside the differentiable geometry; solver norms weight the rows.

Cost. Element blocks one variable larger, sources evaluated at a second set of points, and a residual whose switches had to be smoothed where the model has kinks. One model only.

Context. With an entropy-conservative volume flux the scheme is provably entropy stable only if the interface flux is. The classical Roe flux is not; es_roe, added later, is, and has the same wave-by-wave dissipation.

Decision. Keep roe as the default and offer es_roe. Measured on the development machine before deciding: on four subsonic presets the two fluxes needed the same Newton steps and Krylov iterations to within 4 % and gave outputs within the discretization error; the transonic preset, which converges with roe in 238 of its 250 steps, did not converge within 250 with es_roe; the residual kernel of es_roe took 1.04 times as long as that of roe on a 2-D case at $p = 3$ and 1.67 times on a 3-D case at $p = 2$, where faces carry more of the work (working tree before revision f36e10e, RTX 3060 Ti).

Cost. The default configuration carries no entropy proof, which every page that mentions stability has to say. Revisit when the transonic case converges with the entropy-stable flux.

Context. Total derivatives are exact up to the tolerance of a linear solve. A solve that stops short returns numbers that look like derivatives.

Decision. solve_linear, the total products and totals raise LinearSolveError with the reduction they achieved; returning the result with a flag is something the caller asks for. The OpenMDAO components turn it into an analysis error.

Cost. Cases that used to return a number now stop. The first was the adjoint of the turbulent flat plate with block Jacobi alone (what it showed).

Context. On the wall-resolved mesh of the turbulent flat plate, block Jacobi alone left the adjoint solve short of its tolerance (decision 12 made that visible), and the piecewise-constant multigrid level only repaired it. The standard remedy for such meshes is to solve implicitly along lines of strongly coupled elements. Gyre assembles no matrix, so the usual implementation, which picks the line blocks out of a stored Jacobian, was not available.

Decision. Build the lines from the geometry of the mesh (face area over centre distance, a threshold on the anisotropy, links only where both elements agree), take the two off-diagonal blocks of a line element from derivatives of the split-argument residual with respect to the neighbour state, with one batch of tangents per colour of a colouring of the line elements, and invert the block-tridiagonal systems exactly. The lines do not depend on the flow, the same code serves every multigrid level and the transposed solve, and a mesh without anisotropy gets block Jacobi back unchanged. The option is solver.smoother: line; the default smoother and the preset flat_plate_sa are unchanged.

Cost. Three more batches of block derivatives at every update of the preconditioner (three colours on aligned lines), two more blocks of storage per stored line slot (how the slots are bounded), and a longer first compilation. Two comparisons are recorded, the turbulent flat plate in two dimensions and a laminar sphere in three. In both, the linear solves to a tight tolerance need fewer than half the Krylov iterations (study), and in both the flow solve from the free stream, whose loose Newton steps spend their time in assembling the preconditioner, gains nothing:

Making it the default for wall-resolved cases needs the same measurement on a large mesh and on a turbulent case in three dimensions, which has not been made.

Context. In three dimensions the element blocks of the smoother, and the working memory of their assembly, decide what fits on one device. Everything else the solver keeps is a few dozen state vectors.

Decision. solver.precond_precision: single rounds the blocks, column by column as they are assembled, stores them and applies them in single precision. The inverses and the factorisation of the line systems are still computed in double precision, because that is where a badly conditioned block loses digits, and the operator, the Krylov vectors and the transfers are never touched: a rounded preconditioner changes how fast flexible GMRES converges, not what it converges to. The option is part of the case, so a session, the error estimate and the adjoint use it consistently.

Cost and evidence. On the three cases of the comparison the iteration counts and the answers do not change. The largest 3-D solve is run with it. The default remains double: three small cases and one large one are not every flow, shocks with artificial viscosity among them, and a default that changes the iteration path of every existing case should be changed on more evidence than that.

No optimizer, no design-variable abstraction, no geometry engine: those belong to the MDO framework. No bespoke workflow scheduler beyond the small pipeline engine. No database: the store is files.