---
title: Viscous terms — the BR2 scheme
summary: Viscous fluxes use the second scheme of Bassi and Rebay with lifting operators evaluated on Gauss–Lobatto nodes. The discretization is compact (face neighbours only) and symmetric, which the implicit solver and the adjoint both rely on.
section: theory
order: 5
level: 3
keywords: [BR2, lifting operator, interior penalty, adjoint consistency, viscous flux]
---
<!-- gyre docs: expanded Markdown of theory/viscous-br2.md (generated tables, verification numbers and chart data are written out). Base level of this page: L3 Mechanics; passages at another level are bracketed by 'begin level N' / 'end level N' comments and start with [LN Name]. Built from revision ceb1f4b0be. -->

A discontinuous solution has no gradient at element faces, and using only element-interior
gradients gives an inconsistent — and unstable — diffusion operator. Gyre uses the BR2 scheme:
the gradient is corrected by *lifting* the solution jumps at faces into the elements.

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

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:

$$
\mathbf r_f = \frac{s_J}{w_\text{end}\, J}\,(\mathbf W^* - \mathbf W^-)\otimes\mathbf n,
\qquad
\mathbf W^* = \begin{cases} \tfrac12(\mathbf W^- + \mathbf W^+) & \text{interior face} \\ \mathbf W_b & \text{boundary face.} \end{cases}
$$

The scheme then uses two different gradients:

$$
\underbrace{\mathbf q = \mathbf G + \sum_f \mathbf r_f}_{\text{volume flux } \mathbf F_v(\mathbf U,\mathbf q)}
\qquad\qquad
\underbrace{\mathbf q_f = \mathbf G + \eta\,\mathbf r_f}_{\text{face flux } \tfrac12\left[\mathbf F_v(\mathbf U^-, \mathbf q_f^-) + \mathbf F_v(\mathbf U^+, \mathbf q_f^+)\right]\cdot\mathbf n}
$$

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.

<!-- end level 3 -->

<!-- begin level 4 (Theory) -->
**[L4 Theory] Why this pairing is symmetric**

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

$$
\sum_K (\mathbf K\,\mathbf R([\mathbf W]),\, \mathbf G(\mathbf v))_K = -\sum_f \langle [\mathbf W],\, \{\!\{\mathbf K^{\mathsf T}\mathbf G(\mathbf v)\}\!\}\cdot\mathbf n\rangle_f ,
$$

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](error-estimation.md) 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.

**Table: Relative asymmetry and smallest eigenvalue** (Structure of the BR2 viscous operator)

| warp | boundaries | p | asym. (heat) | asym. (stress) | min eig (heat) | min eig (stress) |
|---|---|---|---|---|---|---|
| 0 | periodic | 2 | 2.7046e-16 | 2.0254e-16 | -5.4470e-16 | 6.5560e-16 |
| 0 | periodic | 3 | 2.2839e-16 | 2.5421e-16 | 1.1123e-16 | -1.8926e-15 |
| 0 | Dirichlet | 2 | 7.9058e-17 | 1.1709e-16 | 0.0423156 | 0.0252442 |
| 0 | Dirichlet | 3 | 6.5254e-17 | 1.4417e-16 | 0.0246366 | 0.0147011 |
| 0.05 | periodic | 2 | 1.9043e-16 | 2.8210e-16 | 1.5375e-15 | 5.5369e-17 |
| 0.05 | periodic | 3 | 2.3864e-16 | 2.9273e-16 | 1.1192e-15 | -5.0847e-16 |
| 0.05 | Dirichlet | 2 | 5.4250e-17 | 1.5964e-16 | 0.0417777 | 0.0247796 |
| 0.05 | Dirichlet | 3 | 1.3399e-16 | 1.6369e-16 | 0.0242548 | 0.0143994 |

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

<!-- end level 4 -->

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

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.

<!-- end level 3 -->

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

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](implicit-solver.md).

**Verification status: Structure of the BR2 viscous operator** (study `viscous_operator`): **PASSED**, 17 of 17 checks passed, 116.9 s.

At a quiescent state the energy–energy and momentum–momentum blocks of the weak-form Jacobian are exactly the discrete heat-conduction and viscous-stress operators. Their symmetry (the hallmark of an adjoint-consistent diffusion discretization) and positive semi-definiteness are checked on Cartesian and warped meshes with periodic and Dirichlet boundaries, and the Jacobian sparsity pattern is checked to be compact (face neighbours only), which element-block preconditioning relies on.

| Check | Value | Expected | Result |
|---|---|---|---|
| warp=0.0, periodic, p=2: asymmetry | 2.7046e-16 | < 1e-12 | pass |
| warp=0.0, periodic, p=2: smallest eigenvalue | -5.4470e-16 | >= -1e-12 | pass |
| warp=0.0, periodic, p=3: asymmetry | 2.5421e-16 | < 1e-12 | pass |
| warp=0.0, periodic, p=3: smallest eigenvalue | -1.8926e-15 | >= -1e-12 | pass |
| warp=0.0, Dirichlet, p=2: asymmetry | 1.1709e-16 | < 1e-12 | pass |
| warp=0.0, Dirichlet, p=2: smallest eigenvalue | 0.0252442 | >= -1e-12 | pass |
| warp=0.0, Dirichlet, p=3: asymmetry | 1.4417e-16 | < 1e-12 | pass |
| warp=0.0, Dirichlet, p=3: smallest eigenvalue | 0.0147011 | >= -1e-12 | pass |
| warp=0.05, periodic, p=2: asymmetry | 2.8210e-16 | < 1e-12 | pass |
| warp=0.05, periodic, p=2: smallest eigenvalue | 5.5369e-17 | >= -1e-12 | pass |
| warp=0.05, periodic, p=3: asymmetry | 2.9273e-16 | < 1e-12 | pass |
| warp=0.05, periodic, p=3: smallest eigenvalue | -5.0847e-16 | >= -1e-12 | pass |
| warp=0.05, Dirichlet, p=2: asymmetry | 1.5964e-16 | < 1e-12 | pass |
| warp=0.05, Dirichlet, p=2: smallest eigenvalue | 0.0247796 | >= -1e-12 | pass |
| warp=0.05, Dirichlet, p=3: asymmetry | 1.6369e-16 | < 1e-12 | pass |
| warp=0.05, Dirichlet, p=3: smallest eigenvalue | 0.0143994 | >= -1e-12 | pass |
| element couplings outside the face-neighbour stencil | 0 | == 0 | pass |

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

<!-- end level 3 -->

<!-- begin level 2 (Practice) -->
**[L2 Practice] What accuracy to expect**

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](../verification/order-of-accuracy.md)). 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.**

**Table: Observed orders on the finest mesh pair** (Manufactured solutions: Navier–Stokes)

| regime | p | L2 error (finest) | state | velocity | pressure | temperature |
|---|---|---|---|---|---|---|
| heat conduction (quiescent) | 2 | 3.0899e-05 | 3.60088 | 3.6433 | 3.60046 | 3.05238 |
| heat conduction (quiescent) | 3 | 9.7099e-07 | 3.99849 | 4.94591 | 3.99795 | 4.58993 |
| diffusion dominated, Re = 10 | 2 | 4.9402e-04 | 2.33256 | 2.72687 | 2.33663 | 3.1517 |
| diffusion dominated, Re = 10 | 3 | 2.7784e-06 | 3.79396 | 4.19354 | 3.79776 | 4.61291 |
| diffusion dominated, no mean flow, Re = 10 | 2 | 0.00135697 | 1.99132 | 2.57181 | 1.99047 | 2.70647 |
| convection dominated, Re = 1000 | 2 | 8.6201e-05 | 2.852 | 2.99477 | 2.84167 | 2.95746 |
| convection dominated, Re = 1000 | 3 | 8.8491e-07 | 3.91862 | 3.98075 | 3.90575 | 4.11096 |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T10:48:07-0500; device gpu · NVIDIA GeForce RTX 3060 Ti; 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: [mms_navier_stokes.json](../data/verification/mms_navier_stokes.json).

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Why the rate is reduced: a hypothesis, with the evidence for and against**

*Hypothesis.* When diffusion dominates, the momentum and energy equations behave like an
equal-order discretization of a Stokes-like system, in which the pressure-like variables are
controlled only through the velocity divergence and lose one order.

*Evidence for.*

* The diffusion operator alone reaches $p+1$: heat conduction in a quiescent gas, and an
  independent one-dimensional BR2 model that includes the nonlinear work term.
* The viscous operator is symmetric to round-off, so the loss is not an adjoint-inconsistency
  of the implementation.
* Pressure is the slowest variable and the loss shrinks as convection takes over, which is the
  pattern the hypothesis predicts.

*Evidence against, and what is missing.*

* In the compressible equations solved here pressure is a thermodynamic variable, not a
  Lagrange multiplier of a divergence constraint. The manufactured solution runs at Mach 0.5,
  far from the incompressible limit in which the Stokes argument is exact.
* Nothing in the repository isolates the mechanism. Three other candidates are present in this
  very discretization and have not been excluded: under-integration of the BR2 form by
  Gauss–Lobatto quadrature, the penalty sitting exactly on the coercivity bound (see above),
  and the non-conservative gradient in $(\mathbf u, T)$.
* The studies that would discriminate — the rate against Mach number, Gauss against
  Gauss–Lobatto quadrature, a sweep of $\eta$ — have not been run.

Until they are, treat the rates in the table as the fact and the explanation as a conjecture.

<!-- end level 3 -->

<!-- end level 2 -->