---
title: Nonlinearly Stable Flux Reconstruction (NSFR)
url: https://www.emergentmind.com/topics/nonlinearly-stable-flux-reconstruction-nsfr
type: topic
---

# Nonlinearly Stable Flux Reconstruction (NSFR)

Nonlinearly stable flux reconstruction (NSFR) is a class of high-order flux-reconstruction discretizations in which the FR correction is embedded into a modified mass matrix and the nonlinear volume terms are written in split form so that the semidiscrete scheme satisfies a discrete entropy inequality. In the compressible-flow setting, NSFR merges the energy-stable FR framework with entropy-stable DG constructions through skew-symmetric stiffness operators, entropy-projected states, and entropy-conservative or entropy-stable two-point fluxes [2103.02647][2312.07725]. Subsequent work extended the same core formulation to curvilinear grids, fully discrete relaxation Runge–Kutta time integration, viscous turbulence and LES/ILES, adaptive lifting for shocks, and space-time discretizations [2409.00208][2411.12108][2511.01564][2604.19900].

## 1. Research lineage and conceptual position

A precursor to later NSFR formulations showed that flux reconstruction can be recast into the residual-distribution framework and vice versa, and used that connection to give a first demonstration of entropy stability for the FR schemes under consideration on general polygonal meshes [1807.01261]. The 2021 split-form development then derived, for the first time, nonlinearly stable ESFR schemes in split form for uncollocated, modal ESFR discretizations with different volume and surface cubature nodes; the key enabling step was applying the splitting to the discrete stiffness operator rather than only to face liftings [2103.02647]. The 2023 curvilinear extension generalized the construction to three-dimensional compressible flow on curvilinear coordinates and introduced a weight-adjusted modified mass matrix, while the 2024 fully discrete work coupled the semidiscretization with relaxation Runge–Kutta to control temporal entropy production [2312.07725][2409.00208].

The 2024 LES study shifted the emphasis from formal nonlinear stability to under-resolved viscous turbulence, showing how the FR correction parameter can be used to stabilize and accelerate implicit LES without over-integration [2411.12108]. Two 2025 branches then addressed shock-dominated regimes: one introduced an adaptive choice of lifting operator based on Persson’s sensor, and the other examined bound-preserving limiters, positivity, and shock-capturing behavior across several FR parameters, quadratures, and two-point fluxes [2511.01564][2507.09131]. In 2026, NSFR was embedded in a fully implicit space-time scheme using FR in space and DG in time, with fully discrete entropy-preservation for \(c_{DG}\) and entropy-stability for small \(c\) [2604.19900].

| arXiv id | Contribution | Setting |
|---|---|---|
| [1807.01261] | FR recast as residual distribution; entropy-stable construction | Polygonal meshes |
| [2103.02647] | First split-form nonlinearly stable ESFR derivation | 1D nonlinear conservation laws |
| [2312.07725] | Curvilinear, weight-adjusted NSFR | 3D compressible flow |
| [2409.00208] | Fully-discrete NSFR via relaxation RK | Burgers, Euler, Navier–Stokes |
| [2411.12108] | LES/ILES assessment of NSFR | Viscous free-shear turbulence |
| [2511.01564] | Adaptive lifting operator | Shock robustness |
| [2507.09131] | Positivity limiter and shock-capturing study | 1D/2D Euler |
| [2604.19900] | Space-time NSFR | Fully implicit space-time FR/DG |

This lineage places NSFR at the intersection of FR, ESFR/VCJH correction theory, SBP split forms, and entropy-stable DG. A plausible implication is that NSFR is best understood not as a separate discretization family unrelated to DG, but as a filtered DG-type formulation whose distinctive feature is the replacement of the standard mass inverse by \((M+K)^{-1}\) together with split-form flux differencing.

## 2. Semidiscrete formulation and the role of the FR correction parameter

For viscous compressible flow, NSFR is formulated for the compressible Navier–Stokes equations in conservation form,
$$
\frac{\partial \mathbf U}{\partial t}
+ \frac{\partial \mathbf F^c_i(\mathbf U)}{\partial x_i}
- \frac{\partial \mathbf F^v_i(\mathbf U,\nabla \mathbf U)}{\partial x_i}
= \mathbf 0,
$$
with \(\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T\), and in the nondimensional form the viscous stress tensor is
$$
\tau_{ij} = \frac{2\mu}{\Re_\infty}
\Bigl(S_{ij}-\tfrac13 S_{kk}\delta_{ij}\Bigr),
\qquad
S_{ij}=\tfrac12\bigl(\partial_j v_i+\partial_i v_j\bigr)
$$
[2411.12108].

The domain is typically partitioned into hexahedral elements. On each element \(m\), the discrete solution and its gradient are expanded in a tensor-product modal basis of order \(p\),
$$
u_m^h(\mathbf x,t)=\sum_{i=1}^{(p+1)^3}\hat u_{m,i}(t)\,\chi_{m,i}(\mathbf x),
\qquad
q_m^h(\mathbf x,t)=\sum_i\hat q_{m,i}(t)\,\chi_{m,i}(\mathbf x),
$$
and the strong-form semidiscretization is written compactly as
$$
(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt}
=
-\bigl[\Vol-\Surf\{\cdots\}\bigr]
-\bigl[\Visc-\Surf\{\cdots\}\bigr].
$$
Here \(\mathbf K_m\) is the FR correction matrix and depends on a scalar parameter \(c\) [2411.12108].

That parameter indexes the FR family. In the LES-oriented and shock-oriented formulations, \(c=c_{DG}\) recovers DG, \(c=c_+\) is the maximum-filter FR, and intermediate values such as \(c_{HU}\) reproduce Huynh-type corrections [2409.00208][2507.09131]. In one dimension the modified matrix can be written as
$$
K = c (D^p)^T M D^p,
$$
while in multidimensions it is assembled from tensor-product derivative operators [2409.00208]. In the curvilinear framework the same idea appears as a modified mass matrix \(M_m+K_m\), with a weight-adjusted approximation of its inverse to avoid dense per-element inversions [2312.07725].

On general polygonal meshes, the older RD-FR formulation expresses the FR correction in Raviart–Thomas spaces \(RT_p(K)\), with correction functions chosen so that their normal traces enforce the face-flux conditions and the residuals remain conservative [1807.01261]. This suggests that NSFR is not restricted conceptually to tensor-product hexahedra, even though the most developed compressible-flow formulations exploit tensor-product structure for efficiency.

## 3. Entropy stability, split forms, and limits of exact invariant preservation

The defining analytical property of NSFR is nonlinear stability in the entropy sense. The volume term is written using a skew-symmetric split form, often through Chan’s hybrid operator \(\widetilde Q-\widetilde Q^T\), and paired with an entropy-conservative two-point flux satisfying Tadmor’s shuffle condition,
$$
\bigl(\mathbf v_L-\mathbf v_R\bigr)^T \mathbf f^*(u_L,u_R)
=
\psi(u_L)-\psi(u_R).
$$
Under these conditions, the semidiscrete scheme satisfies a discrete entropy inequality of the form
$$
\sum_m d_t\!\int \eta(\mathbf u_m^h)\,d\Omega \le 0.
$$
A central result quoted in the viscous LES study states: “For any \(c\ge c_{\min}\), the FR scheme with split-form volume terms and an EC two-point flux is provably entropy-stable without added dissipation or de-aliasing” [2411.12108].

The 2021 split-form derivation clarified why earlier ESFR split formulations could fail nonlinearly. In NSFR, the ESFR filter \((M+K)^{-1}\) is applied to both conservative and non-conservative volume integrals and to face liftings, so that the SBP identity is realized at the fully discrete algebraic level for uncollocated modal or nodal bases [2103.02647]. The 2023 curvilinear paper extended that logic to mapped elements satisfying the discrete geometric conservation law, and also established exact free-stream preservation and global conservation under the discrete GCL [2312.07725].

A recurrent misconception is that entropy stability implies machine-precision kinetic-energy preservation on arbitrary quadratures. The curvilinear analysis proves the opposite: “A high-order scheme cannot be discretely kinetic-energy preserving to machine precision whenever surface quadrature \(\not\subset\) volume quadrature” [2312.07725]. In practice, collocated LGL settings preserve the discrete KE balance to \(O(10^{-12})\) in the reported tests, whereas GL settings drift [2312.07725].

The viscous turbulence study exposed a related diagnostic issue. For the Taylor–Green vortex, the pressure-dilatation-based dissipation rate is consistent with literature when obtained from the kinetic-energy budget terms, but direct computation of \(\langle P\,(\nabla\!\cdot v)\rangle\) exhibits spurious oscillations. Those oscillations are significantly lower for a collocated scheme and are effectively eliminated by adding Roe upwind dissipation to the two-point numerical flux, leading the authors to attribute them to the treatment of face terms in nonlinearly stable schemes [2411.12108]. This does not contradict entropy stability; rather, it distinguishes entropy bounds from the pointwise smoothness of individual post-processed budget terms.

## 4. Implementation, quadrature, and computational structure

The practical efficiency of NSFR relies on tensor-product implementation. In the LES formulation, all tensor-product operations—mass-matrix multiplies, differentiation matrices, and Hadamard products in the split-form volume term—are evaluated by sum-factorization, reducing the per-element cost from \(O(p^6)\) to \(O(p^4)\) [2411.12108]. The curvilinear formulation states the same idea in dimension-agnostic form: a three-dimensional \(n^3\times n^3\) matrix-vector product is replaced by three one-dimensional \(n\times n\) products at cost \(O(d\,n^{d+1})\) rather than \(O(n^{2d})\), with correspondingly lower storage [2312.07725].

For curved elements, the modified mass inverse \((M_m+K_m)^{-1}\) becomes dense. The weight-adjusted form approximates it through reference-element operators and elementwise multiplication by \((WJ_m)^{-1}\), so that no dense element matrices are formed [2312.07725]. This is one of the main enablers for extending NSFR beyond affine meshes without losing the low-storage character associated with FR and DGSEM implementations.

Quadrature choice affects both robustness and invariant properties. The 2023 curvilinear study verified exact discrete entropy balance to \(O(10^{-13})\) for both GL and LGL quadratures, but only the collocated LGL case achieved machine-precision KE conservation in the reported TGV test [2312.07725]. In the 2025 shock study, GLL quadrature was preferred because it allowed larger CFL values and simpler face corrections, whereas GL produced slightly sharper smooth solutions but forced CFL \(\approx 0.05\) in shocked tests to avoid entropy-projection spikes [2507.09131].

Relative to classical DG with over-integration, NSFR was consistently reported as cheaper in the tested regimes. In the viscous turbulence comparison on 16 MPI ranks, Figure 21 showed about a \(2\times\) speedup of NSFR relative to DG with over-integration for \(p\in[5,12]\) [2411.12108]. In the 3D curvilinear TGV study on eight AMD cores, overintegrated DG was \(3\)–\(4\times\) more expensive than NSFR-EC, while NSFR-EC cost only about \(10\)–\(15\%\) extra CPU time relative to standard conservative DG, yet conservative DG diverged at \(p=5\) whereas NSFR-EC remained robust [2312.07725].

## 5. Turbulence, LES, and implicit LES behavior

The most detailed assessment of NSFR for viscous turbulence is the study of subsonic free-shear flow and the Taylor–Green vortex. For DNS verification at \(\Re=1600\), a \(p=7\), \(32^3\)-element NSFR.IR-GL discretization with \(c=c_{DG}\) and \(256^3\) DOF reproduced kinetic energy, dissipation, enstrophy, and \(\varepsilon_P\) in line with a \(512^3\) finite-difference reference [2411.12108]. Under-resolved runs then showed stable and accurate implicit LES without explicit SGS modeling or over-integration: \(p=5\) at \(96^3\) DOF, and \(p=7\) at \(64^3\) DOF, both remained stable, with the latter capturing the transitional phase well [2411.12108].

A central numerical observation was that increasing the FR correction parameter enlarges explicit time-step limits while preserving the desired ILES behavior. With \(c=c_{DG}\), NSFR.IR-GL used \(\mathrm{CFL}\simeq 0.26\), compared with \(\mathrm{CFL}\simeq 0.14\) for DG with over-integration; pushing to \(c=c_+\) increased the limit to \(\mathrm{CFL}\simeq 0.36\), with all runs advanced by SSP-RK3 [2411.12108]. The same study reported that the choice of entropy-conservative two-point flux—IR, KG, CH, or \(CH_{RA}\)—produced virtually identical TGV/LES results, so the two-point flux choice did not affect stability for that problem [2411.12108]. That conclusion is therefore problem-specific, not a universal statement about shock-dominated flows.

The kinetic-energy budget was analyzed through
$$
\frac{dE_K}{dt}
=
\underbrace{\langle P\,\partial_i v_i\rangle}_{+\varepsilon_P}
-
\underbrace{\frac{2}{\Re_\infty}\langle \mu\,S^d_{ij}S^d_{ij}\rangle}_{+\varepsilon_v}
\equiv \varepsilon_P-\varepsilon_v,
$$
and the “observed” \(\varepsilon_P\) was also reconstructed from \(\varepsilon-\varepsilon_v\) in post-processing [2411.12108]. Standard Smagorinsky, shear-improved Smagorinsky, dynamic Smagorinsky, and their high-pass-filtered versions did not improve ILES accuracy for TGV; they added dissipative bias at small scales [2411.12108]. This is a direct challenge to the assumption that explicit eddy-viscosity closures are automatically beneficial once the spatial discretization is entropy stable.

Spectral diagnostics required additional care. Accurate TKE spectra demanded oversampling of the continuous velocity field onto an equi-spaced grid of \(2(p+1)\) points per direction inside each element before FFTs; without oversampling, a spurious TKE pile-up appeared near the cutoff wavenumber \(\kappa_c\) [2411.12108]. For decaying homogeneous isotropic turbulence in the Comte-Bellot–Corrsin configuration at \(\Re_\lambda\sim 70\), an NSFR \(p=3\), \(128^3\)-DOF discretization captured the experimental spectra at \(t^*=1,2\) with less than \(2\%\) error in total \(E_K\) [2411.12108].

## 6. Shock capturing, positivity preservation, and adaptive lifting

In shock-dominated Euler problems, NSFR is typically paired with a positivity or bound-preserving limiter. One 2025 study extended the Zhang–Shu limiter by including the minimum density and pressure at the solution nodes when computing the rescaling parameter, so that positivity is monitored not only on mixed tensor quadrature grids but also at the nodal set itself [2507.09131]. The same work emphasized that NSFR combines provable nonlinear stability with the increased time step from energy-stable FR, and found that stronger FR filtering improves robustness, CFL limits, and mitigation of overshoots and oscillations [2507.09131].

The reported parameter study distinguished sharply between turbulence and shock regimes. In shocked tests, the \(CH_{RA}\) two-point flux was found most robust, KG under-dissipated and could violate entropy stability, and GLL quadrature was preferred in most cases [2507.09131]. The one-dimensional FR parameter family was reported as \(c_{DG}=0\), \(c_{SD}\approx 7.4\times 10^{-6}\), \(c_{HU}\approx 1.3\times 10^{-5}\), \(c_+\approx 2.9\times 10^{-5}\), and \(c_{+\times 10}=2.9\times 10^{-4}\) [2507.09131]. In the Sod tube at \(p=3\), \(N=512\), NSFR with \(c_{DG}\) survived without limiter up to \(\mathrm{CFL}=0.2\), and with the positivity-preserving limiter ran at \(\mathrm{CFL}=0.5\), while doubling \(c\to c_+\) reduced overshoot near the shock by about \(30\%\) [2507.09131]. In the Leblanc tube, \(c_+\) ran at \(\mathrm{CFL}=0.3\) without TVD; in 2D vortex–shock interaction, NSFR with \(c_{DG}\) was stable to \(\mathrm{CFL}=0.27\) without limiter and \(\mathrm{CFL}=0.58\) with limiter [2507.09131].

A complementary 2025 note made the FR correction parameter adaptive. It introduced Persson’s modal sensor
$$
S_e=\frac{\sum_{|\ell|=p}\hat u_\ell^2}{\sum_{|\ell|\le p}\hat u_\ell^2},
\qquad
s_e=\log_{10}S_e,
\qquad
s_0=-4\log_{10}p,
$$
followed by a smooth cutoff and the local choice
$$
c_{\rm elem}=\varepsilon\,c_+.
$$
The method reverts toward DG in smooth regions and turns on FR-level dissipation near shocks [2511.01564]. On a Gaussian-pulse convergence problem, DG and FR both converged at order \(p+1\), FR had larger \(L^1\), \(L^2\), and \(L^\infty\) error constants than DG, and the adaptive scheme lay strictly between them while recovering full \(p+1\) order [2511.01564]. On the Leblanc shock tube with 1920 DOF, the maximum stable CFLs were \(0.10\), \(0.005\), and failure for DG at \(p=3,4,5\); \(0.29\), \(0.21\), and \(0.15\) for the adaptive scheme; and \(0.30\), \(0.23\), and \(0.19\) for FR [2511.01564]. In 2D shock diffraction at \(p=3\), the maximum CFL was \(0.45\) for DG and \(0.61\) for both adaptive and FR [2511.01564].

These shock results qualify a common simplification about NSFR. The scheme by itself does not eliminate oscillations in the way a dedicated shock-capturing method does; the adaptive-lifting paper states explicitly that it cannot eliminate such oscillations, but together with a positivity-preserving limiter it provides solutions that are essentially oscillation-free [2511.01564]. The practical shock-capturing behavior therefore arises from the coupling of entropy-stable NSFR, FR-parameter tuning, and limiter design.

## 7. Fully discrete and space-time extensions

The 2024 fully-discrete formulation extended entropy-stable NSFR in time using relaxation Runge–Kutta. The update is written as
$$
U^{n+1}=U^n+\gamma^n\Delta t\sum_i b_i f(U^{(i)}),
$$
where \(\gamma^n\) is chosen either from an algebraic formula in the inner-product entropy case or by solving a scalar root-finding problem for a general convex entropy [2409.00208]. For inner-product entropies, FD-NSFR prevents temporal numerical entropy change in the broken Sobolev norm; for general convex numerical entropies it prevents temporal numerical entropy change in the physical \(L_2\) norm, and fully-discrete entropy stability in \(L_2\) is obtained only with the DG correction \(c=0\) because the \(K\)-contribution otherwise leaves an \(O(\Delta t^p)\) remainder [2409.00208].

The numerical consequences were problem dependent. For inviscid Burgers, FD-NSFR preserved energy to machine zero while the semidiscrete entropy-stable method displayed \(O(\Delta t^p)\) drift, and \(\gamma\to 1\) at \(O(\Delta t^{p-1})\) [2409.00208]. For inviscid Taylor–Green vortex, the reported maximum stable CFL for the semidiscrete scheme was \(0.48\) for \(c=0\) and \(0.54\) for \(c_+\); FD-NSFR showed zero numerical entropy growth to machine precision, while the semidiscrete method showed \(O(10^{-7})\)–\(O(10^{-8})\) drift [2409.00208]. The same study reported that FD-NSFR required about one eighth as many time steps to satisfy the same “entropy-conserving” criterion, with a per-step overhead of about \(1.3\)–\(1.6\times\) relative to the semidiscrete scheme [2409.00208].

The 2026 space-time formulation replaced method-of-lines by FR in space and DG in time on tensor-product space-time elements. In that framework, \(c=0\) recovers space-time DGSEM, \(c=c_{Hu}\) recovers Huynh-type FR, and \(c\to\infty\) approaches the spectral-difference method [2604.19900]. The space-time nonlinearly stable FR scheme uses skew-symmetric stiffness operators in both space and time and is fully-discretely entropy preserving with \(c_{DG}\) or entropy-stable for small \(c\) [2604.19900]. Numerically, linear advection and Euler tests confirmed \(p+1\) convergence for \(c\le c_{Hu}\), with a drop to order \(p\) once \(c\) exceeds a threshold near \(c_{Hu}\); the abstract reports a reduction in computational cost up to about \(70\%\) as \(c\) is increased [2604.19900].

Taken together, these developments indicate that NSFR has evolved from a semidiscrete entropy-stable reinterpretation of ESFR into a broader framework covering curvilinear geometry, viscous turbulence, shock-adaptive filtering, fully discrete temporal stabilization, and fully implicit space-time discretization. The persistent structural theme is the same in every variant: a filtered mass operator \(M+K\), split-form flux differencing, and entropy-compatible interelement coupling.

Source: https://www.emergentmind.com/topics/nonlinearly-stable-flux-reconstruction-nsfr