Nonlinearly Stable Flux Reconstruction (NSFR)

• NSFR is a high-order discretization method that embeds flux reconstruction corrections into a modified mass matrix and employs split forms to enforce discrete entropy inequalities.
• It leverages a tunable correction parameter to interpolate between DG and Huynh-type schemes, balancing numerical dissipation and stability in compressible flow simulations.
• The framework extends to curvilinear grids, fully discrete time integration, turbulence modeling, and shock capturing through adaptive lifting and bound-preserving limiters.

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 (Cicchino et al., 2021, Cicchino et al., 2023). 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 (Pethrick et al., 2024, Brillon et al., 2024, Srinivasan et al., 3 Nov 2025, Pethrick et al., 21 Apr 2026).

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 (Abgrall et al., 2018). 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 (Cicchino et al., 2021). 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 (Cicchino et al., 2023, Pethrick et al., 2024).

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 (Brillon et al., 2024). 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 (Srinivasan et al., 3 Nov 2025, Srinivasan et al., 12 Jul 2025). 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$ (Pethrick et al., 21 Apr 2026).

arXiv id	Contribution	Setting
(Abgrall et al., 2018)	FR recast as residual distribution; entropy-stable construction	Polygonal meshes
(Cicchino et al., 2021)	First split-form nonlinearly stable ESFR derivation	1D nonlinear conservation laws
(Cicchino et al., 2023)	Curvilinear, weight-adjusted NSFR	3D compressible flow
(Pethrick et al., 2024)	Fully-discrete NSFR via relaxation RK	Burgers, Euler, Navier–Stokes
(Brillon et al., 2024)	LES/ILES assessment of NSFR	Viscous free-shear turbulence
(Srinivasan et al., 3 Nov 2025)	Adaptive lifting operator	Shock robustness
(Srinivasan et al., 12 Jul 2025)	Positivity limiter and shock-capturing study	1D/2D Euler
(Pethrick et al., 21 Apr 2026)	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)$$

(Brillon et al., 2024).

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 $c$0 is the FR correction matrix and depends on a scalar parameter $c$1 (Brillon et al., 2024).

That parameter indexes the FR family. In the LES-oriented and shock-oriented formulations, $c$2 recovers DG, $c$3 is the maximum-filter FR, and intermediate values such as $c$4 reproduce Huynh-type corrections (Pethrick et al., 2024, Srinivasan et al., 12 Jul 2025). In one dimension the modified matrix can be written as

$c$5

while in multidimensions it is assembled from tensor-product derivative operators (Pethrick et al., 2024). In the curvilinear framework the same idea appears as a modified mass matrix $c$6, with a weight-adjusted approximation of its inverse to avoid dense per-element inversions (Cicchino et al., 2023).

On general polygonal meshes, the older RD-FR formulation expresses the FR correction in Raviart–Thomas spaces $c$7, with correction functions chosen so that their normal traces enforce the face-flux conditions and the residuals remain conservative (Abgrall et al., 2018). 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 $c$8, and paired with an entropy-conservative two-point flux satisfying Tadmor’s shuffle condition,

$c$9

Under these conditions, the semidiscrete scheme satisfies a discrete entropy inequality of the form

$(M+K)^{-1}$0

A central result quoted in the viscous LES study states: “For any $(M+K)^{-1}$1, the FR scheme with split-form volume terms and an EC two-point flux is provably entropy-stable without added dissipation or de-aliasing” (Brillon et al., 2024).

The 2021 split-form derivation clarified why earlier ESFR split formulations could fail nonlinearly. In NSFR, the ESFR filter $(M+K)^{-1}$2 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 (Cicchino et al., 2021). 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 (Cicchino et al., 2023).

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 $(M+K)^{-1}$3 volume quadrature” (Cicchino et al., 2023). In practice, collocated LGL settings preserve the discrete KE balance to $(M+K)^{-1}$4 in the reported tests, whereas GL settings drift (Cicchino et al., 2023).

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 $(M+K)^{-1}$5 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 (Brillon et al., 2024). 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 $(M+K)^{-1}$6 to $(M+K)^{-1}$7 (Brillon et al., 2024). The curvilinear formulation states the same idea in dimension-agnostic form: a three-dimensional $(M+K)^{-1}$8 matrix-vector product is replaced by three one-dimensional $(M+K)^{-1}$9 products at cost $$\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,$$0 rather than $$\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,$$1, with correspondingly lower storage (Cicchino et al., 2023).

For curved elements, the modified mass inverse $$\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,$$2 becomes dense. The weight-adjusted form approximates it through reference-element operators and elementwise multiplication by $$\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,$$3, so that no dense element matrices are formed (Cicchino et al., 2023). 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 $$\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,$$4 for both GL and LGL quadratures, but only the collocated LGL case achieved machine-precision KE conservation in the reported TGV test (Cicchino et al., 2023). 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 $$\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,$$5 in shocked tests to avoid entropy-projection spikes (Srinivasan et al., 12 Jul 2025).

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 1 showed about a $$\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,$$6 speedup of NSFR relative to DG with over-integration for $$\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,$$7 (Brillon et al., 2024). In the 3D curvilinear TGV study on eight AMD cores, overintegrated DG was $$\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,$$8–$$\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,$$9 more expensive than NSFR-EC, while NSFR-EC cost only about $\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T$0–$\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T$1 extra CPU time relative to standard conservative DG, yet conservative DG diverged at $\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T$2 whereas NSFR-EC remained robust (Cicchino et al., 2023).

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 $\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T$3, a $\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T$4, $\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T$5-element NSFR.IR-GL discretization with $\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T$6 and $\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T$7 DOF reproduced kinetic energy, dissipation, enstrophy, and $\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T$8 in line with a $\mathbf U=[\rho,\rho v_1,\rho v_2,\rho v_3,E]^T$9 finite-difference reference (Brillon et al., 2024). Under-resolved runs then showed stable and accurate implicit LES without explicit SGS modeling or over-integration: $$\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)$$0 at $$\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)$$1 DOF, and $$\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)$$2 at $$\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)$$3 DOF, both remained stable, with the latter capturing the transitional phase well (Brillon et al., 2024).

A central numerical observation was that increasing the FR correction parameter enlarges explicit time-step limits while preserving the desired ILES behavior. With $$\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)$$4, NSFR.IR-GL used $$\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)$$5, compared with $$\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)$$6 for DG with over-integration; pushing to $$\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)$$7 increased the limit to $$\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)$$8, with all runs advanced by SSP-RK3 (Brillon et al., 2024). The same study reported that the choice of entropy-conservative two-point flux—IR, KG, CH, or $$\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)$$9—produced virtually identical TGV/LES results, so the two-point flux choice did not affect stability for that problem (Brillon et al., 2024). That conclusion is therefore problem-specific, not a universal statement about shock-dominated flows.

The kinetic-energy budget was analyzed through

$m$0

and the “observed” $m$1 was also reconstructed from $m$2 in post-processing (Brillon et al., 2024). 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 (Brillon et al., 2024). 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 $m$3 points per direction inside each element before FFTs; without oversampling, a spurious TKE pile-up appeared near the cutoff wavenumber $m$4 (Brillon et al., 2024). For decaying homogeneous isotropic turbulence in the Comte-Bellot–Corrsin configuration at $m$5, an NSFR $m$6, $m$7-DOF discretization captured the experimental spectra at $m$8 with less than $m$9 error in total $p$0 (Brillon et al., 2024).

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 (Srinivasan et al., 12 Jul 2025). 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 (Srinivasan et al., 12 Jul 2025).

The reported parameter study distinguished sharply between turbulence and shock regimes. In shocked tests, the $p$1 two-point flux was found most robust, KG under-dissipated and could violate entropy stability, and GLL quadrature was preferred in most cases (Srinivasan et al., 12 Jul 2025). The one-dimensional FR parameter family was reported as $p$2, $p$3, $p$4, $p$5, and $p$6 (Srinivasan et al., 12 Jul 2025). In the Sod tube at $p$7, $p$8, NSFR with $p$9 survived without limiter up to $$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),$$0, and with the positivity-preserving limiter ran at $$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),$$1, while doubling $$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),$$2 reduced overshoot near the shock by about $$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),$$3 (Srinivasan et al., 12 Jul 2025). In the Leblanc tube, $$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),$$4 ran at $$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),$$5 without TVD; in 2D vortex–shock interaction, NSFR with $$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),$$6 was stable to $$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),$$7 without limiter and $$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),$$8 with limiter (Srinivasan et al., 12 Jul 2025).

A complementary 2025 note made the FR correction parameter adaptive. It introduced Persson’s modal sensor

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

followed by a smooth cutoff and the local choice

$(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt} = -\bigl[\Vol-\Surf\{\cdots\}\bigr] -\bigl[\Visc-\Surf\{\cdots\}\bigr].$0

The method reverts toward DG in smooth regions and turns on FR-level dissipation near shocks (Srinivasan et al., 3 Nov 2025). On a Gaussian-pulse convergence problem, DG and FR both converged at order $(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt} = -\bigl[\Vol-\Surf\{\cdots\}\bigr] -\bigl[\Visc-\Surf\{\cdots\}\bigr].$1, FR had larger $(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt} = -\bigl[\Vol-\Surf\{\cdots\}\bigr] -\bigl[\Visc-\Surf\{\cdots\}\bigr].$2, $(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt} = -\bigl[\Vol-\Surf\{\cdots\}\bigr] -\bigl[\Visc-\Surf\{\cdots\}\bigr].$3, and $(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt} = -\bigl[\Vol-\Surf\{\cdots\}\bigr] -\bigl[\Visc-\Surf\{\cdots\}\bigr].$4 error constants than DG, and the adaptive scheme lay strictly between them while recovering full $(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt} = -\bigl[\Vol-\Surf\{\cdots\}\bigr] -\bigl[\Visc-\Surf\{\cdots\}\bigr].$5 order (Srinivasan et al., 3 Nov 2025). On the Leblanc shock tube with 1920 DOF, the maximum stable CFLs were $(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt} = -\bigl[\Vol-\Surf\{\cdots\}\bigr] -\bigl[\Visc-\Surf\{\cdots\}\bigr].$6, $(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt} = -\bigl[\Vol-\Surf\{\cdots\}\bigr] -\bigl[\Visc-\Surf\{\cdots\}\bigr].$7, and failure for DG at $(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt} = -\bigl[\Vol-\Surf\{\cdots\}\bigr] -\bigl[\Visc-\Surf\{\cdots\}\bigr].$8; $(\mathbf M_m+\mathbf K_m)\,\frac{d\hat{\mathbf u}_m}{dt} = -\bigl[\Vol-\Surf\{\cdots\}\bigr] -\bigl[\Visc-\Surf\{\cdots\}\bigr].$9, $c$00, and $c$01 for the adaptive scheme; and $c$02, $c$03, and $c$04 for FR (Srinivasan et al., 3 Nov 2025). In 2D shock diffraction at $c$05, the maximum CFL was $c$06 for DG and $c$07 for both adaptive and FR (Srinivasan et al., 3 Nov 2025).

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 (Srinivasan et al., 3 Nov 2025). 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

$c$08

where $c$09 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 (Pethrick et al., 2024). 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 $c$10 norm, and fully-discrete entropy stability in $c$11 is obtained only with the DG correction $c$12 because the $c$13-contribution otherwise leaves an $c$14 remainder (Pethrick et al., 2024).

The numerical consequences were problem dependent. For inviscid Burgers, FD-NSFR preserved energy to machine zero while the semidiscrete entropy-stable method displayed $c$15 drift, and $c$16 at $c$17 (Pethrick et al., 2024). For inviscid Taylor–Green vortex, the reported maximum stable CFL for the semidiscrete scheme was $c$18 for $c$19 and $c$20 for $c$21; FD-NSFR showed zero numerical entropy growth to machine precision, while the semidiscrete method showed $c$22–$c$23 drift (Pethrick et al., 2024). 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 $c$24–$c$25 relative to the semidiscrete scheme (Pethrick et al., 2024).

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$26 recovers space-time DGSEM, $c$27 recovers Huynh-type FR, and $c$28 approaches the spectral-difference method (Pethrick et al., 21 Apr 2026). 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$29 or entropy-stable for small $c$30 (Pethrick et al., 21 Apr 2026). Numerically, linear advection and Euler tests confirmed $c$31 convergence for $c$32, with a drop to order $c$33 once $c$34 exceeds a threshold near $c$35; the abstract reports a reduction in computational cost up to about $c$36 as $c$37 is increased (Pethrick et al., 21 Apr 2026).

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 $c$38, split-form flux differencing, and entropy-compatible interelement coupling.