Viscous terms — the BR2 scheme

Let $\mathbf W = (\mathbf u, T)$ be the gradient variables and $\mathbf G = \tfrac1J \sum_a (J\mathbf a^a)\,\mathbf D_a \mathbf W$ the element-local (strong, non-conservative) gradient. For each face $f$ of an element define the face lifting, non-zero only at that face’s nodes:

The scheme then uses two different gradients:

and the divergence of the volume flux is taken in weak form with $\mathbf D^w = \mathbf M^{-1} \mathbf D^{\mathsf T}\mathbf M$, the adjoint of $\mathbf D$ in the quadrature inner product. Using only the face’s own lifting in the face flux (rather than the full $\mathbf q$, as BR1 does) is what keeps the stencil compact.

Because the LGL mass matrix is diagonal, the lifting relation $\int_K \mathbf r_f\cdot\boldsymbol\tau = \int_f (\mathbf W^* - \mathbf W^-)\,\mathbf n\cdot\boldsymbol\tau$ is satisfied pointwise by the formula above — no linear solve. Summing the volume term over the two elements adjacent to a face gives, for a linear flux $\mathbf F_v = \mathbf K\mathbf q$,

which is the mirror image of the consistency term $-\sum_f\langle\{\!\{\mathbf K\mathbf G(\mathbf W)\}\!\}\cdot\mathbf n, [\mathbf v]\rangle_f$ produced by the face flux, and the penalty term $\eta\sum_f\langle\mathbf r_f([\mathbf W]),\mathbf r_f([\mathbf v])\rangle$ is symmetric positive semi-definite. The discrete operator is therefore the symmetric BR2 bilinear form evaluated with LGL quadrature. Symmetry is not cosmetic: it is the discrete statement of adjoint consistency, the property that lets functionals converge faster than the solution and makes the output error estimate sharp.

The requirement that the gradient ($\mathbf D$, strong, non-conservative) and divergence ($\mathbf D^w$, weak) be exact adjoints is why Gyre mixes the two forms rather than using one throughout.

The coercivity proofs for BR2 (Brezzi, Manzini, Marini, Pietra and Russo 2000; Arnold, Brezzi, Cockburn and Marini 2002) need $\eta$ strictly greater than the number of faces of an element; at equality the argument gives only a non-negative, not a coercive, form. Gyre’s default is the customary borderline value $\eta = 2d$ (4 for quadrilaterals, 6 for hexahedra), which lies on the edge of the proven range, not inside it. The proofs are also for exact integration, whereas Gyre uses Gauss–Lobatto quadrature. No instability has been observed with the default in the tests and studies; the value is an input of the discretization (Discretization(..., eta_br2=...) in Python) for anyone who wants to be inside the proven range. At the finite-volume level $p = 0$ there is no interior gradient and the lifting alone carries the viscous flux. With $\eta = 2$ it is the two-point difference along the face normal: exact on uniform spacing, too large by $(h_1 + h_2)^2 / (4 h_1 h_2)$ between cells of widths $h_1$ and $h_2$. The gradient components tangential to a face, and with them the cross-derivative terms of the stress, are absent. The $p = 0$ viscous operator is therefore a p-multigrid preconditioner level, not a consistent Navier–Stokes discretization, and cannot be selected as the solution degree of a viscous case.

The residual of an element depends only on itself and its face neighbours. This is verified by inspecting the sparsity of a dense AD Jacobian, and it is what allows the element-diagonal blocks used by the preconditioner to be extracted exactly with the split-argument trick.

The diffusion operator on its own reaches the design rate: heat conduction in a quiescent gas converges at $p+1$, and an independent one-dimensional implementation of exactly this formulation converges at $p+1$ for every degree (study). For the coupled compressible Navier–Stokes system in a diffusion-dominated regime the state converges at a reduced rate between $p$ and $p+1$, with pressure the slowest variable. Convection-dominated flows recover approximately $p+1$. The reduced rate is a measurement; its cause has not been established.