---
title: Implicit solver
summary: Steady states are found by pseudo-transient Newton–Krylov iteration. The Jacobian is never formed — its action comes from automatic differentiation — and it is preconditioned by p-multigrid with element-block-Jacobi smoothing whose blocks are also obtained by AD.
section: theory
order: 7
level: 2
keywords: [Newton-Krylov, pseudo-transient continuation, GMRES, p-multigrid, block Jacobi, line smoother, line-implicit, matrix-free, mixed precision, single precision]
---
<!-- gyre docs: expanded Markdown of theory/implicit-solver.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. -->

The steady problem is $\mathbf R_w(\mathbf U) = \mathbf 0$, where $\mathbf R_w = (wJ)\,\mathbf R$ is the
weak-form residual (the pointwise residual scaled by the diagonal mass matrix). Scaling by the
mass balances large and small elements in every norm the solver uses, and makes adjoint
variables nodal values of the continuous adjoint.

## Pseudo-transient continuation

Each nonlinear iteration solves

$$
\left( \frac{\mathbf M}{\Delta\tau} + \frac{\partial \mathbf R_w}{\partial \mathbf U} \right) \Delta\mathbf U = -\mathbf R_w(\mathbf U),
\qquad \Delta\tau_e = \mathrm{CFL}\cdot\Delta t_e ,
$$

with an element-local time scale $\Delta t_e$ from the convective and viscous spectral radii.
Small CFL gives robust, damped steps far from the solution; as CFL grows the method becomes
Newton's method with an exact Jacobian. The linear systems are solved to a fixed relative
tolerance (inexact Newton, `lin_tol`, $10^{-2}$ by default) and the CFL grows by at most a
factor 10 per step, so the terminal convergence is linear, not quadratic: the residual falls
by about the factor `lin_tol` per step. The
[solver study](../verification/solvers.md) checks that contraction.

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Globalisation**

* **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](../guide/python-api.md)).

**Chart: Residual history** (Nonlinear and linear solver convergence). History chart of |Rw| against nonlinear iteration, y axis log; 6 series. The table holds the plotted data.

| series | nonlinear iteration | \|Rw\| |
|---|---|---|
| Euler NACA 0012, p = 2: block Jacobi | 0 | 0.366317 |
| Euler NACA 0012, p = 2: block Jacobi | 1 | 0.238266 |
| Euler NACA 0012, p = 2: block Jacobi | 2 | 0.153701 |
| Euler NACA 0012, p = 2: block Jacobi | 3 | 0.128176 |
| Euler NACA 0012, p = 2: block Jacobi | 4 | 0.107526 |
| Euler NACA 0012, p = 2: block Jacobi | 5 | 0.0911261 |
| Euler NACA 0012, p = 2: block Jacobi | 6 | 0.0708977 |
| Euler NACA 0012, p = 2: block Jacobi | 7 | 0.0507337 |
| Euler NACA 0012, p = 2: block Jacobi | 8 | 0.0305722 |
| Euler NACA 0012, p = 2: block Jacobi | 9 | 0.0140312 |
| Euler NACA 0012, p = 2: block Jacobi | 10 | 0.00539924 |
| Euler NACA 0012, p = 2: block Jacobi | 11 | 9.9688e-04 |
| Euler NACA 0012, p = 2: block Jacobi | 12 | 6.0154e-05 |
| Euler NACA 0012, p = 2: block Jacobi | 13 | 8.0257e-07 |
| Euler NACA 0012, p = 2: block Jacobi | 14 | 7.8549e-09 |
| Euler NACA 0012, p = 2: block Jacobi | 15 | 7.7948e-11 |
| Euler NACA 0012, p = 2: block Jacobi | 16 | 1.0461e-12 |
| Euler NACA 0012, p = 2: p-multigrid | 0 | 0.366317 |
| Euler NACA 0012, p = 2: p-multigrid | 1 | 0.238446 |
| Euler NACA 0012, p = 2: p-multigrid | 2 | 0.15374 |
| Euler NACA 0012, p = 2: p-multigrid | 3 | 0.127815 |
| Euler NACA 0012, p = 2: p-multigrid | 4 | 0.10725 |
| Euler NACA 0012, p = 2: p-multigrid | 5 | 0.0910524 |
| Euler NACA 0012, p = 2: p-multigrid | 6 | 0.0707134 |
| Euler NACA 0012, p = 2: p-multigrid | 7 | 0.0507297 |
| Euler NACA 0012, p = 2: p-multigrid | 8 | 0.0306123 |
| Euler NACA 0012, p = 2: p-multigrid | 9 | 0.0141012 |
| Euler NACA 0012, p = 2: p-multigrid | 10 | 0.00544363 |
| Euler NACA 0012, p = 2: p-multigrid | 11 | 0.00100033 |
| Euler NACA 0012, p = 2: p-multigrid | 12 | 6.0277e-05 |
| Euler NACA 0012, p = 2: p-multigrid | 13 | 6.6459e-07 |
| Euler NACA 0012, p = 2: p-multigrid | 14 | 6.4517e-09 |
| Euler NACA 0012, p = 2: p-multigrid | 15 | 4.5689e-11 |
| Euler NACA 0012, p = 2: p-multigrid | 16 | 7.0588e-13 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 0 | 20.2116 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 1 | 13.1423 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 2 | 9.47857 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 3 | 6.10869 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 4 | 4.09279 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 5 | 3.38947 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 6 | 1.97246 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 7 | 0.483745 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 8 | 0.141892 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 9 | 0.074575 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 10 | 0.0199166 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 11 | 0.0024487 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 12 | 3.2684e-05 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 13 | 3.2557e-07 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 14 | 3.2121e-09 |
| Laminar cylinder Re = 40, p = 2: block Jacobi | 15 | 3.0809e-11 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 0 | 20.2116 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 1 | 13.1323 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 2 | 9.44805 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 3 | 6.06778 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 4 | 4.06751 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 5 | 3.34489 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 6 | 1.94072 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 7 | 0.477141 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 8 | 0.141383 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 9 | 0.0735677 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 10 | 0.0196538 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 11 | 0.00245037 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 12 | 2.8469e-05 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 13 | 1.9200e-07 |
| Laminar cylinder Re = 40, p = 2: p-multigrid | 14 | 1.6586e-09 |
| Euler NACA 0012, p = 3: block Jacobi | 0 | 0.311716 |
| Euler NACA 0012, p = 3: block Jacobi | 1 | 0.313325 |
| Euler NACA 0012, p = 3: block Jacobi | 2 | 0.138704 |
| Euler NACA 0012, p = 3: block Jacobi | 3 | 0.0970839 |
| Euler NACA 0012, p = 3: block Jacobi | 4 | 0.074958 |
| Euler NACA 0012, p = 3: block Jacobi | 5 | 0.0617388 |
| Euler NACA 0012, p = 3: block Jacobi | 6 | 0.0472716 |
| Euler NACA 0012, p = 3: block Jacobi | 7 | 0.0336438 |
| Euler NACA 0012, p = 3: block Jacobi | 8 | 0.0197131 |
| Euler NACA 0012, p = 3: block Jacobi | 9 | 0.0086172 |
| Euler NACA 0012, p = 3: block Jacobi | 10 | 0.00352864 |
| Euler NACA 0012, p = 3: block Jacobi | 11 | 6.8180e-04 |
| Euler NACA 0012, p = 3: block Jacobi | 12 | 4.6864e-05 |
| Euler NACA 0012, p = 3: block Jacobi | 13 | 6.2728e-07 |
| Euler NACA 0012, p = 3: block Jacobi | 14 | 6.1512e-09 |
| Euler NACA 0012, p = 3: block Jacobi | 15 | 6.0011e-11 |
| Euler NACA 0012, p = 3: block Jacobi | 16 | 1.0292e-12 |
| Euler NACA 0012, p = 3: p-multigrid | 0 | 0.311716 |
| Euler NACA 0012, p = 3: p-multigrid | 1 | 0.313727 |
| Euler NACA 0012, p = 3: p-multigrid | 2 | 0.13878 |
| Euler NACA 0012, p = 3: p-multigrid | 3 | 0.0967875 |
| Euler NACA 0012, p = 3: p-multigrid | 4 | 0.074779 |
| Euler NACA 0012, p = 3: p-multigrid | 5 | 0.0616541 |
| Euler NACA 0012, p = 3: p-multigrid | 6 | 0.0471596 |
| Euler NACA 0012, p = 3: p-multigrid | 7 | 0.0336347 |
| Euler NACA 0012, p = 3: p-multigrid | 8 | 0.0197123 |
| Euler NACA 0012, p = 3: p-multigrid | 9 | 0.00862647 |
| Euler NACA 0012, p = 3: p-multigrid | 10 | 0.00352607 |
| Euler NACA 0012, p = 3: p-multigrid | 11 | 6.7471e-04 |
| Euler NACA 0012, p = 3: p-multigrid | 12 | 4.6216e-05 |
| Euler NACA 0012, p = 3: p-multigrid | 13 | 5.4278e-07 |
| Euler NACA 0012, p = 3: p-multigrid | 14 | 3.3127e-09 |
| Euler NACA 0012, p = 3: p-multigrid | 15 | 1.8404e-11 |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T10:42:49-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: [solver_convergence.json](../data/verification/solver_convergence.json).

<!-- end level 3 -->

## Matrix-free Krylov solver

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Linearize once, apply many times**

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*:

```python
_, 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.

<!-- end level 3 -->

## Preconditioning without assembling a matrix

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Exact element blocks from split arguments**

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:

$$
\mathbf B_e = \frac{\mathbf M_e}{\Delta\tau_e} + \left.\frac{\partial \mathbf R_{w,e}}{\partial \mathbf U_e}\right|_{\text{neighbours frozen}} .
$$

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

<!-- end level 3 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] p-multigrid**

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.

**Table: Solver effort** (Nonlinear and linear solver convergence)

| case | preconditioner | levels | Newton steps | Krylov iterations | max Krylov / step | converged | final / reference residual | terminal contraction (mean) | terminal contraction (worst step) | wall time [s] |
|---|---|---|---|---|---|---|---|---|---|---|
| Euler NACA 0012, p = 2 | block Jacobi | 2 | 16 | 1206 | 182 | yes | 2.8557e-12 | 0.0109235 | 0.0134203 | 58.5489 |
| Euler NACA 0012, p = 2 | p-multigrid | 2, 1, 0 | 16 | 142 | 17 | yes | 1.9270e-12 | 0.010203 | 0.0154498 | 104.262 |
| Laminar cylinder Re = 40, p = 2 | block Jacobi | 2 | 15 | 1707 | 320 | yes | 1.5243e-12 | 0.00980497 | 0.00996111 | 48.199 |
| Laminar cylinder Re = 40, p = 2 | p-multigrid | 2, 1, 0 | 14 | 143 | 19 | yes | 8.2062e-11 | 0.0076328 | 0.00863853 | 73.4588 |
| Euler NACA 0012, p = 3 | block Jacobi | 3 | 16 | 1360 | 202 | yes | 3.3017e-12 | 0.0117945 | 0.0171502 | 74.6616 |
| Euler NACA 0012, p = 3 | p-multigrid | 3, 1, 0 | 15 | 154 | 23 | yes | 5.9042e-11 | 0.00582304 | 0.00610324 | 79.4693 |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T10:42:49-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: [solver_convergence.json](../data/verification/solver_convergence.json).

<!-- end level 3 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Line smoothing on anisotropic meshes**

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):

$$
\mathbf T = \begin{pmatrix}
\mathbf B_1 & \mathbf U_1 \\
\mathbf L_2 & \mathbf B_2 & \mathbf U_2 \\
 & \ddots & \ddots & \ddots \\
 & & \mathbf L_m & \mathbf B_m
\end{pmatrix},
\qquad
\mathbf L_j = \frac{\partial \mathbf R_{w,e_j}}{\partial \mathbf U_{e_{j-1}}},\quad
\mathbf U_j = \frac{\partial \mathbf R_{w,e_j}}{\partial \mathbf U_{e_{j+1}}} .
$$

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

**Table: Linear solves at the converged flow (tolerance 1e-10)** (Line-implicit smoothing on a wall-resolved mesh)

| preconditioner | system | Krylov iterations | residual reduction | tolerance reached | wall time of the Krylov solve, compiled, preconditioner assembled before [s] |
|---|---|---|---|---|---|
| block Jacobi, levels 2-1-0 | forward | 238 | 8.9804e-11 | yes | 6.18439 |
| block Jacobi, levels 2-1-0 | transposed (adjoint) | 247 | 9.4097e-11 | yes | 10.9862 |
| line Jacobi, levels 2-1-0 | forward | 66 | 9.5595e-11 | yes | 4.98718 |
| line Jacobi, levels 2-1-0 | transposed (adjoint) | 70 | 9.4827e-11 | yes | 7.34169 |
| block Jacobi, one level | forward | 720 | 1.1248e-04 | no | 4.00668 |
| block Jacobi, one level | transposed (adjoint) | 720 | 2.8702e-04 | no | 12.8319 |
| line Jacobi, one level | forward | 441 | 9.9505e-11 | yes | 5.62206 |
| line Jacobi, one level | transposed (adjoint) | 465 | 8.7870e-11 | yes | 5.46822 |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T10:58:43-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: [line_smoother.json](../data/verification/line_smoother.json).

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](../verification/solvers.md#line-smoothing-on-a-wall-resolved-mesh)
records the wall times of two small cases, one in two and one in three dimensions; nothing
is claimed beyond them.

**Table: One assembly of the preconditioner at the converged flow** (Line-implicit smoothing on a wall-resolved mesh)

| preconditioner | wall time, compiled, fastest of five [s] |
|---|---|
| block Jacobi, levels 2-1-0 | 0.1873 |
| line Jacobi, levels 2-1-0 | 4.04077 |
| block Jacobi, one level | 3.00486 |
| line Jacobi, one level | 4.08774 |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T10:58:43-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: [line_smoother.json](../data/verification/line_smoother.json).

<!-- end level 3 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Single-precision blocks**

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:

**Table: Double against single precision of the smoother's blocks** (Single-precision blocks in the preconditioner)

| case | blocks | unknowns | nonlinear steps | Krylov iterations, flow solve | drag coefficient | Krylov iterations, forward solve | Krylov iterations, transposed solve | linear solutions: relative difference to double | stored blocks [MB] |
|---|---|---|---|---|---|---|---|---|---|
| turbulent flat plate (SA) | double | 6750 | 35 | 540 | 0.00575365 | 238 | 249 | — | 2.94 |
| turbulent flat plate (SA) | single | 6750 | 35 | 540 | 0.00575365 | 238 | 249 | 1.7075e-11 | 1.47 |
| turbulent flat plate (SA), line smoother | double | 6750 | 35 | 267 | 0.00575365 | 66 | 63 | — | 8.8984 |
| turbulent flat plate (SA), line smoother | single | 6750 | 35 | 267 | 0.00575365 | 66 | 63 | 4.4650e-11 | 4.4492 |
| laminar sphere (3-D), Re = 100, degree 1 | double | 3840 | 14 | 79 | 2.13486 | 46 | 43 | — | 1.248 |
| laminar sphere (3-D), Re = 100, degree 1 | single | 3840 | 14 | 79 | 2.13486 | 46 | 43 | 4.1034e-14 | 0.624 |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T11:27:02-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: [mixed_precision.json](../data/verification/mixed_precision.json).

What it buys is in the record of the
[largest solve](../verification/three-dimensional.md#largest-solve-on-one-device), which
states what happened when the same problem was tried with double-precision blocks on the
same device. The default stays `double` ([decision 14](../engineering/design-decisions.md)).

<!-- end level 3 -->

<!-- begin level 4 (Theory) -->
**[L4 Theory] Line smoothing in detail: lines, colours, factors**

**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,

$$
\Bigl(\frac{\partial \mathbf R_{w}}{\partial \mathbf U_\text{nbr}}\,\mathbf t_c\Bigr)_e
  = \sum_{f \in \mathcal N(e),\ \operatorname{colour}(f) = c}
    \frac{\partial \mathbf R_{w,e}}{\partial \mathbf U_f}\,\mathbf v .
$$

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

$$
\mathbf c = \mathbf D^{-1}\mathbf b,\qquad
\mathbf z_j = \mathbf c_j - \mathbf G_j \mathbf z_{j-1},\qquad
\mathbf x_j = \mathbf z_j - \mathbf W_j \mathbf x_{j+1},
$$

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

$$
\mathbf y_j = \mathbf b_j - \mathbf W_{j-1}\T \mathbf y_{j-1},\qquad
\mathbf u_j = \mathbf y_j - \mathbf G_{j+1}\T \mathbf u_{j+1},\qquad
\mathbf x = \mathbf D^{-\mathsf T}\mathbf u .
$$

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.

<!-- end level 4 -->

<!-- begin level 4 (Theory) -->
**[L4 Theory] Cost model and its consequence**

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](../verification/solvers.md#what-the-iterations-cost)). 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](../verification/solvers.md) records Krylov iteration counts and does not
show that the V-cycle lowers total work or wall time.

<!-- end level 4 -->