Implicit solver

• CFL evolution. After an accepted full step the CFL is multiplied by the residual reduction factor, clipped to $[1.5, 10]$; after a rejected step it is divided by 10.
• Line search. A step length $\omega \in \{1, \tfrac12, \dots\}$ is accepted if the state stays admissible ($\rho, p > 0$) and the unsteady residual $\|\mathbf R_w(\mathbf U + \omega\Delta\mathbf U) + \tfrac{\mathbf M}{\Delta\tau}\omega\Delta\mathbf U\|$ drops below $\|\mathbf R_w(\mathbf U)\|$. Requiring the steady residual to drop would reject legitimate transient steps at low CFL.
• Inexact Newton. Linear systems are solved to a relative tolerance of $10^{-2}$ by default.
• Warm starts. When the initial state is the converged solution of a neighbouring problem — the normal situation inside a design loop — continuation starts directly at $\mathrm{CFL} = 10^4$; the line search still protects against a bad step.
• Round-off floor. If the residual stops halving for six steps in the Newton regime, the solver stops and reports whether the plateau is below the acceptance level.
• Failure is reported, not raised. A solve ends unconverged when the step limit is reached, the residual stalls above the acceptance level, the residual is not finite or the pseudo-time step collapses; the history carries converged: false and the reason. The discipline then keeps its last converged state as the warm start of the next solve, holds the failed iterate separately (failed_state), and refuses derivative calls until a solve converges (Python API).

The linear solver is restarted flexible GMRES written in JAX control flow (jax.lax.while_loop). The operator is built inside the same traced function:

_, lin = jax.linearize(lambda V: residual_weak(V), U)   # one primal sweep, kept by XLA
A = lambda v: idm * v + lin(v)                          # tangent sweep per Krylov iteration
x, info = fgmres(A, M, rhs, x0, m=40, restarts=2, tol=1e-2)

Because linearisation and Krylov loop live in one compiled graph, the primal intermediate values are computed once and every iteration costs only a tangent evaluation. The adjoint solver is the same code with jax.vjp in place of jax.linearize: the transposed operator is exact, with no separately derived adjoint residual.

Orthogonalisation is classical Gram–Schmidt with one re-orthogonalisation pass — two dense matrix–vector products per iteration, as robust as modified Gram–Schmidt and far better suited to GPUs.

The residual kernel takes its state twice, residual(U_own, geom, th, U_nbr): U_own supplies everything element-local, U_nbr supplies every neighbour trace (and neighbour gradients). Evaluated at U_own = U_nbr = U it is the residual; differentiated with respect to U_own alone it is exactly block diagonal, one dense block per element:

A single JVP whose tangent is the same unit vector in every element therefore returns one column of all blocks simultaneously. All blocks cost block_size batched JVPs ($n_\text{var}(p+1)^d$, e.g. 64 for $p=3$ in 2-D), independent of the number of elements, with no graph colouring and no sparse assembly. The blocks are inverted in one batched call. The linearization study multiplies these inverses with the diagonal blocks of a dense AD Jacobian and finds the identity to $10^{-13}$.

Block Jacobi removes high-frequency, element-local error but not smooth error, so its Krylov iteration counts grow with mesh size. Gyre wraps it in a V-cycle over polynomial degrees (default $p, \lfloor p/2\rfloor, \dots, 0$) on the same mesh:

• levels: full discretizations at lower degree, linearised about the $L^2$-projected state; the lowest level $p=0$ is a finite-volume scheme on the curved cells;
• smoother: damped block Jacobi ($\omega = 0.7$, two pre- and post-sweeps) or a few block-Jacobi-preconditioned GMRES iterations;
• transfers: prolongation is nodal interpolation $\mathbf P$; residual restriction is $\mathbf P^{\mathsf T}$ — the Galerkin restriction, exact because the residual is in weak form;
• adjoint: the transposed hierarchy uses the same transfers with transposed level operators and transposed blocks, so forward and adjoint solves share one implementation.

A wall-resolved mesh has cells that are hundreds or thousands of times longer than they are thick. Across the thin direction an element is coupled to its neighbours far more strongly than in the others, and a smoother that freezes the neighbours, as block Jacobi does, leaves the error along that direction almost untouched. The coarser levels of the hierarchy have the same cells and inherit the anisotropy: the piecewise-constant level makes the linear solves on such a mesh converge, but they stay slow (first rows of the table below).

With smoother: line the elements are grouped into lines that follow the strong coupling, and the smoother inverts the block-tridiagonal system of every line exactly (Mavriplis 1998; Fidkowski, Oliver, Lu and Darmofal 2005):

• Lines come from the mesh. The coupling of two elements that share a face is the face area divided by the distance of their centres, the coefficient of a two-point diffusion flux. An element joins a line through the pair of opposite faces of its strongest direction if that direction is at least four times stronger than the next one (for a box-shaped cell: an aspect ratio of two), and two elements are linked only where both choose the shared face. Elements of isotropic regions belong to no line and are smoothed by block Jacobi; on a mesh without lines the option changes nothing. The lines do not depend on the flow and are the same on every level.
• Off-diagonal blocks from the neighbour argument. The derivative of the split-argument residual with respect to U_nbr contains every coupling block. The line elements are coloured so that a tangent which is one unit vector in all elements of one colour puts into each line element a column of the block of exactly one line neighbour. Aligned lines need three colours, so the two off-diagonal blocks cost three times the diagonal ones.
• Factorisation. $\mathbf D_1 = \mathbf B_1$, $\mathbf D_j = \mathbf B_j - \mathbf L_j \mathbf D_{j-1}^{-1} \mathbf U_{j-1}$, followed by one forward and one backward substitution per application. All lines are processed together; only the steps along a line are sequential.
• Adjoint. The transposed systems are solved with the same factors, $\mathbf T\T = (\mathbf I + \mathbf W)\T (\mathbf D + \mathbf L)\T$ with $\mathbf W_j = \mathbf D_j^{-1}\mathbf U_j$, so forward and adjoint solves again share one implementation. The unit tests compare both with a dense solve of $\mathbf T$ cut out of an AD Jacobian, in two and three dimensions and with the turbulence model.

What the option costs. The assembly of the preconditioner takes one more batch of block derivatives per colour on every level, three on aligned lines, so about four times the derivatives of block Jacobi; the measured time of one assembly is in the table below. The storage is the block-Jacobi inverse of every element plus two blocks per stored line slot; the lines are stored in a few tables by length, so the slots exceed the line elements only by the padding within a table (the summary of a case and gyre mesh info report both numbers). The study records the wall times of two small cases, one in two and one in three dimensions; nothing is claimed beyond them.

The blocks of the smoother are the largest thing the solver stores: $n_\text{var}^2 (p+1)^{2d}$ numbers per element and level, and two more blocks per line slot with the line smoother. In three dimensions they, and the working memory of their assembly, are what limits the size of a problem on one device. They are also the one part of the solver that need not be exact: flexible GMRES converges with any preconditioner, and converges to the solution of the system it is given.

With solver.precond_precision: single

• every column of the blocks is rounded to single precision as the batched tangents produce it, so the blocks never exist in double precision and the assembly needs half the memory as well as the storage;
• the inverses, and the recurrence that factorises the line systems, are computed in double precision from the rounded blocks, a batch of elements at a time, and their results are rounded: a block of a wall-resolved cell with a turbulence model can have a condition number that single precision does not resolve, so the inversion is where the digits are needed;
• the smoother is applied in single precision, to a rounded copy of the vector it is given.

The residual, the Jacobian–vector and vector–Jacobian products, the Krylov basis and the multigrid transfers stay in double precision. The discrete equations, and therefore the converged flow, the adjoint and every derivative, are unchanged; only the number of Krylov iterations could differ, and on the cases measured it does not:

What it buys is in the record of the largest solve, which states what happened when the same problem was tried with double-precision blocks on the same device. The default stays double (decision 14).

Lines. For a face shared by the elements $e$ and $f$ the coupling is $c_{ef} = A_{ef} / (d_e + d_f)$, with $A_{ef}$ the area of the face and $d_e$, $d_f$ the distances of the two element centres from the face centre (so a periodic neighbour is at its true distance); centres and areas are those of the straight-sided element through the corner nodes. A quadrilateral or hexahedron has one pair of opposite faces per local direction $a$; its strength is $s_a(e) = \max(c_{e,2a}, c_{e,2a+1})$. The element chooses its strongest direction if $s_{a^*} \ge 4\, s_{\text{second}}$. For a box of sizes $h_a$, $c \propto h_b h_c / h_a$, so the ratio of two directions is $(h_b/h_a)^2$ and the threshold is an aspect ratio of two. Two elements are linked where each chooses the face they share; every element then has at most two links, so the connected sets are paths, or closed loops across a periodic direction, which are cut at their weakest link. A single anisotropic layer next to isotropic cells forms no line, because its neighbour does not choose the shared face.

Colours. Let $\mathcal N(e)$ be $e$ and its face neighbours that belong to lines. A tangent that is the unit vector $\mathbf v$ in every element of colour $c$ and zero elsewhere gives, differentiated through the neighbour argument,

This is one column of the block of a line neighbour $f$ of $e$ exactly when no other member of $\mathcal N(e)$ shares the colour of $f$: that is the condition the colouring enforces, for every line element and both of its line neighbours. On a block of aligned lines the position along the line modulo three satisfies it (the line neighbours of $e$ have the two other colours, the neighbours across the lines have the colour of $e$ itself), hence three colours. Where lines of different directions meet (the wall-normal lines and the circumferential loops of an O-mesh) more can be needed, and no single heuristic is best on every mesh: a saturation-degree colouring and three greedy passes in different orders are tried, and the colouring with the fewest colours is kept. The number is part of the case summary. The cost of assembling the preconditioner is therefore $(1 + n_\text{colours})\,n_\text{block}$ batched tangents per level where block Jacobi takes $n_\text{block}$.

Factors. Per line, with $\mathbf D_j$ from the recurrence above, $\mathbf G_j = \mathbf D_j^{-1}\mathbf L_j$ and $\mathbf W_j = \mathbf D_j^{-1}\mathbf U_j$ are stored, and $\mathbf D_j^{-1}$ takes the place of the block-Jacobi inverse of the element: two blocks per line slot in addition to the one per element that block Jacobi stores. Lines of very different length occur together (the circumferential loops in the far field of an O-mesh next to its wall-normal lines), and one table as long as the longest line would be mostly padding that is factorised, stored and swept like real slots. The lines are therefore sorted by length and cut into at most four tables so that the stored slots are fewest; a further table is added only while the slots exceed the line elements by more than a quarter. One application is

and the transposed one, from $\mathbf T\T = (\mathbf I + \mathbf W)\T(\mathbf D + \mathbf L)\T$,

The first (last) product is taken for all elements at once and is all that elements outside the lines receive, which is block Jacobi. The two sweeps run, table by table, over the slots of the longest line of the table; each step is one batched block product over its lines.

What stays approximate. $\mathbf T$ contains the coupling along the lines only. The coupling across lines and between elements outside them is left to the Krylov iteration and to the other levels, so line Jacobi is a smoother with the same damping as block Jacobi, not a solver. The lines follow the mesh, not the flow: strong convection across thin cells is not what they are built for. A cell that is thin in two directions at once (the band of cells along the edge where two walls meet in three dimensions) has two strong directions of similar strength, fails the anisotropy test and gets no line: it stays with block Jacobi.

Let $c_R$ be the cost of one residual evaluation. A Krylov iteration at the fine level costs about one tangent sweep ($\approx 1$–$2\,c_R$) plus the V-cycle (a few fine-level operator applications for smoothing plus cheaper coarse levels): counted from the cycle, about 7.6 fine-level applications per iteration with the default sweeps and the levels 2, 1, 0 in two dimensions, against one for block Jacobi. The hierarchy therefore saves operator work only where it cuts the iteration count by more than that factor (measured per case). Building the block inverses costs block_size tangent sweeps of the split residual and is amortised over all Krylov iterations of a Newton step (or several, with precond_lag). Memory is dominated by the Krylov basis (2 × krylov_dim state vectors for flexible GMRES) and the block inverses ($K \cdot$ block_size$^2$): for $p=3$ in 3-D the blocks are $320\times320$ per element, which bounds what fits on a single small GPU.

The dominant fixed cost in practice is JIT compilation: tens of seconds per new combination of mesh size and degree. Sessions exist to pay this once per design study. This cost model is an operation count, not a measurement: the solver study records Krylov iteration counts and does not show that the V-cycle lowers total work or wall time.