discontinuous Galerkin methods, hyperbolic conservation laws, total variation diminishing, total variation bounded, Harten’s lemma, non-SSP time discretizations

1 Introduction↩︎

Discontinuous Galerkin (DG) methods, typically coupled with Runge–Kutta (RK) time discretizations, are widely used for solving hyperbolic conservation laws [1]–[5]. To suppress spurious oscillations near shocks generated by the high-order approximations, total variation diminishing (TVD) or total variation bounded (TVB) limiting is usually applied as a standard remedy [2]. In the classical Runge–Kutta discontinuous Galerkin (RKDG) framework, this limiting strategy is intrinsically tied to what were historically termed TVD RK time discretizations [6], [7], more commonly known today as strong-stability-preserving (SSP) RK time discretizations [8], [9]. These methods can be represented as convex combinations of forward-Euler steps. By ensuring that the forward-Euler step is TVD/TVB in the means, the SSP-RK method inherits this stability property via its underlying convex decomposition. Consequently, the TVD or SSP time integrators were historically believed to be essential for preserving fully discrete TVD properties; indeed, numerical examples in [7] demonstrated that “non-TVD but linearly stable RK time discretization can generate oscillations even for TVD spatial discretization."

To overcome this rigid reliance on SSP structures, we develop a novel TVD/TVB framework tailored for general non-SSP time stepping. Our motivation is twofold. First, explicit SSP-RK methods impose strict order barriers: they require additional stages to achieve fourth-order accuracy and cannot exceed order four with positive weights [7], [10]. Second, many advanced fully discrete DG formulations naturally lack an SSP structure; for instance, the RKDG method with compact stencils (cRKDG) fundamentally requires non-SSP-RK schemes to ensure conservation [11]–[14], while the Arbitrary DERivative (ADER) DG [15]–[18] method does not use RK or multi-step integrators at all, let alone their SSP versions.

Our approach relies on the same theoretical foundation as the classical SSP-RKDG framework, namely Harten’s lemma [19], [20], but applies it from a different perspective. Drawing on this lemma, our key observation is that the TVD/TVB-in-the-means property is dictated solely by the cell-average update in the final corrector stage of the RK integrator. Consequently, this property can be attained by modifying the numerical fluxes in the final stage alone. The interior predictor stages are not constrained by this property; they may use generic approximations without an SSP structure. We leverage this insight by introducing a trace-limited corrector framework for DG schemes paired with generic non-SSP time discretizations. The predictor stages are left unchanged, while their traces entering the final numerical fluxes in the corrector stage are limited with respect to the initial-stage cell averages, thereby ensuring that the final-stage DG solution is TVD/TVB in the means. A standard polynomial limiter is then applied as a postprocessing step to produce the accepted DG polynomial beyond its means. This construction decouples the TVD/TVB limiting mechanism from the SSP structure of the time integrator, while preserving conservation and the usual TVD/TVB-in-the-means interpretation of DG limiting.

As a main application of the proposed framework, we construct the TVD/TVB cRKDG method. The cRKDG method [11] utilizes cell-local derivative operators during the inner stages of an RK method, yielding compact stencils, reduced communication overhead in parallel computing, and simple boundary treatments. However, due to its fundamental reliance on a non-SSP-RK structure, no TVD/TVB cRKDG schemes have been designed in the literature prior to this work. Beyond cRKDG, the proposed framework is also compatible with other fully discrete DG formulations; for instance, Section 3.5 outlines how it can be combined with ADER-DG methods to achieve provable TVD/TVB properties. Moreover, the new framework allows implicit predictors, thereby opening new pathways for designing implicit TVD/TVB DG schemes. To keep our work more focused, we will concentrate on explicit RKDG and cRKDG solvers in subsequent sections; detailed studies on extensions to other fully discrete DG methods will be postponed to future work.

This work aligns with active ongoing efforts in the literature to address two central challenges in numerical conservation laws: oscillation control [21]–[24] and structure preservation [25]–[29]. In particular, we note that the proposed method shares similar goals and outcomes with recent invariant-domain-preserving RK methods (see, e.g., [30]), as both preserve physical invariant properties beyond the SSP structural restriction. However, our approach adopts a fundamentally different methodology. By building upon the classical Harten’s lemma, this approach eliminates the need for immediate limiting operations following each intermediate predictor stage, thereby retaining the compactness of the cRKDG method, streamlining the algorithmic framework for concurrent computation, and facilitating better adaptation to other fully discrete DG formulations beyond the RK family.

The remainder of the paper is organized as follows. Section 2 reviews the classical SSP-RKDG framework and the TVD/TVB polynomial limiter. Section 3 introduces the proposed trace-limited corrector framework for generic RKDG and cRKDG methods, along with a brief discussion on its application to the ADER-DG method. The main properties of the proposed schemes are analyzed in Section 4. Extensions to systems and multidimensional problems are described in Section 5. Numerical tests are presented in Section 6. Conclusions are given in Section 7.

2 Background↩︎

This section reviews the classical SSP-RKDG framework, in which TVD or TVB stability is inherited from the limited forward-Euler DG step.

2.1 Forward-Euler DG methods↩︎

Consider a scalar conservation law in 1D with periodic or compactly supported boundary conditions: \[u_t + f(u)_x = 0, \qquad x\in \Omega, \quad t>0. \label{eq:scl}\tag{1}\]

Let \(\mathcal{T}_h = \{I_i\}_{i=1}^N\), with \(I_i=[x_{i-1/2},x_{i+1/2}]\), be a partition of \(\Omega\). We denote by \(h_i=x_{i+1/2}-x_{i-1/2}\) the cell size of \(I_i\) and \(h = \max_i h_i\). The DG finite element space of degree \(k\) is \(\mathbb{V}_h^k = \bigl\{v\in L^2(\Omega): v|_{I_i}\in \mathbb{P}^k(I_i),\; i=1,\ldots,N\bigr\},\) where \(\mathbb{P}^k(I_i)\) denotes the space of polynomials of degree at most \(k\) on \(I_i\). For \(v\in \mathbb{V}_h^k\), we denote its traces at cell interfaces by \(v_{i+1/2}^\pm = \lim_{\delta\to 0^{+}} v(x_{i+1/2}\pm\delta)\), and write \(\overline{v}_{i} = {1}/{h_i}\cdot\int_{I_i}v(x)\,\mathrm{d}x\) for its cell average on \(I_i\). We also use the shorthands \(\Delta_+\overline{v}_{i} = \overline{v}_{i+1}-\overline{v}_{i}\) and \(\Delta_-\overline{v}_{i} = \overline{v}_{i}-\overline{v}_{i-1}\) for the forward and backward differences of cell averages, respectively.

The standard DG semi-discretization is to find \(u_h(\cdot,t)\in \mathbb{V}_h^k\) such that \[(u_h)_t = \mathcal{F}(u_h), \label{eq:semi-dg}\tag{2}\] where the DG spatial operator \(\mathcal{F}:\mathbb{V}_h^k\to \mathbb{V}_h^k\) is defined by \[\begin{align} \int_{I_i}\mathcal{F}(u_h)v\,\mathrm{d}x = \int_{I_i} f(u_h)v_x\,\mathrm{d}x - \left[ \widehat f(u_{h,i+\frac{1}{2}}^-,u_{h,i+\frac{1}{2}}^+)v_{i+\frac{1}{2}}^- - \widehat f(u_{h,i-\frac{1}{2}}^-,u_{h,i-\frac{1}{2}}^+)v_{i-\frac{1}{2}}^+ \right] \end{align} \label{eq:F-def}\tag{3}\] for all \(v\in \mathbb{V}_h^k\). Here, \(\widehat f\) is a consistent monotone numerical flux, which means \(\widehat f(u,u)=f(u)\), and \(\widehat f\) is nondecreasing in its first argument and nonincreasing in its second argument. Moreover, we assume that \(\widehat f\) is Lipschitz continuous in both arguments.

The raw forward-Euler DG step associated with 2 is \[u_h^{n+1} = u_h^n+\Delta t\,\mathcal{F}(u_h^n). \label{eq:fe-dg}\tag{4}\] Taking \(v=1\) in 3 gives the cell-average update \[\overline{u}_{h,i}^{n+1} = \overline{u}_{h,i}^n -\lambda_i \left[ \widehat f(u_{h,i+1/2}^{n,-},u_{h,i+1/2}^{n,+}) - \widehat f(u_{h,i-1/2}^{n,-},u_{h,i-1/2}^{n,+}) \right], \quad \lambda_i=\frac{\Delta t}{h_i}. \label{eq:fe-cell-average}\tag{5}\] Thus the cell averages evolve by a conservative finite-volume-type update, but with numerical fluxes evaluated at DG traces.

A relevant stability notion is defined through the total variation of cell averages.

The following lemma provides a sufficient condition for TVDM.

2.2 TVD/TVB limiting↩︎

The raw forward-Euler update 4 is not automatically TVDM/TVBM for general DG polynomials, as the interface traces may differ from the cell averages by amounts that violate the Harten form. The classical remedy is to apply a limiter to process the candidate DG solution [2], [19], [31], [32].

The standard minmod function is \[m(a,b,c) = \begin{cases} s\min\{|a|,|b|,|c|\}, &\text{if } \operatorname{sgn}(a)=\operatorname{sgn}(b)=\operatorname{sgn}(c)=s,\\[2mm] 0, &\text{otherwise}. \end{cases} \label{eq:minmod}\tag{6}\] The modified minmod function introduces a parameter \(M\ge 0\) that relaxes 6 by \[\widetilde{m}(a,b,c;M,h) = \begin{cases} a, &\text{if } |a|\le Mh^2,\\[1mm] m(a,b,c), &\text{otherwise}. \end{cases} \label{eq:tvb-minmod}\tag{7}\] The choice \(M=0\) gives the TVD limiting, while \(M>0\) gives the TVB limiting.

We now define the limiter \(\Lambda:\mathbb{V}_h^k\to \mathbb{V}_h^k\). Let \(w\in \mathbb{V}_h^k\) be a candidate DG function. On cell \(I_i\), define the limited traces by \[\tag{8} \begin{align} \widetilde{w}_{i+1/2}^- &= \overline{w}_{i} + \widetilde{m} \left( w_{i+1/2}^- - \overline{w}_{i}, \Delta_+\overline{w}_{i}, \Delta_-\overline{w}_{i}; M,h_i \right), \tag{9} \\ \widetilde{w}_{i-1/2}^+ &= \overline{w}_{i} - \widetilde{m} \left( \overline{w}_{i} - w_{i-1/2}^+, \Delta_+\overline{w}_{i}, \Delta_-\overline{w}_{i}; M,h_i \right). \tag{10} \end{align}\] A cell \(I_i\) is called troubled if at least one of the two evaluations of \(\widetilde{m}\) in 8 modifies its first argument. On cells that are not troubled, the limiter leaves the polynomial unchanged, i.e., \(\Lambda(w)|_{I_i}=w|_{I_i}.\) On a troubled cell \(I_i\), the limiter replaces \(w|_{I_i}\) by the unique reconstructed polynomial in \(\mathbb{P}^{\min(k,2)}(I_i)\) that preserves the cell average while agreeing with the limited trace values in 8 , i.e., \[\overline{\Lambda(w)}_{i} = \overline{w}_i, \qquad \bigl(\Lambda(w)\bigr)_{i+1/2}^{-} = \widetilde{w}_{i+1/2}^{-}, \qquad \bigl(\Lambda(w)\bigr)_{i-1/2}^{+} = \widetilde{w}_{i-1/2}^{+}.\]

Using the convention above, the limited forward-Euler DG method is written as \[w_h^{n+1} = u_h^n+\Delta t\,\mathcal{F}(u_h^n), \qquad u_h^{n+1} = \Lambda(w_h^{n+1}). \label{eq:limited-fe-dg}\tag{11}\] Here \(u_h^n\) is the accepted DG solution at time \(t^n\), while \(w_h^{n+1}\) is the unlimited candidate produced by the forward-Euler step. Since \(\Lambda\) preserves cell averages, the cell-average update of \(u_h^{n+1}\) is the same as that of \(w_h^{n+1}\).

The role of the limiter is to ensure that the traces of the accepted input \(u_h^n\) are compatible with the neighboring cell averages. For the TVD case \(M=0\), under the standard CFL condition associated with the monotone flux, the cell-average update can be written in Harten form ?? –?? . Consequently, the solution to 11 satisfies ?? and is TVDM.

2.3 TVD/TVB stability in SSP-RKDG methods↩︎

SSP-RK methods are designed as convex combinations of forward-Euler steps. Therefore, if the forward-Euler building block is stable in certain seminorms under a time-step restriction, then the SSP method inherits the same stability under the corresponding CFL restriction.

In the RKDG setting, this principle is combined with the limited forward-Euler step 11 for the TVD/TVB stability. For example, the commonly used third-order SSP-RKDG method can be written in the following form: \[\begin{align} w_h^{(1)} &= u_h^n+\Delta t\,\mathcal{F}(u_h^n), & u_h^{(1)} &= \Lambda(w_h^{(1)}), \\ w_h^{(2)} &= u_h^{(1)}+\Delta t\,\mathcal{F}(u_h^{(1)}), & u_h^{(2)} &= \frac{3}{4} u_h^n+ \frac{1}{4}\Lambda(w_h^{(2)}), \\ w_h^{(3)} &= u_h^{(2)}+\Delta t\,\mathcal{F}(u_h^{(2)}), & u_h^{n+1} &= \frac{1}{3} u_h^n+ \frac{2}{3}\Lambda(w_h^{(3)}). \end{align} \label{eq:ssp-rk3-dg}\tag{12}\] Each \(w_h^{(j)}\) is an unlimited forward-Euler candidate, and each \(\Lambda(w_h^{(j)})\) is the corresponding limited DG solution. The accepted stage values are convex combinations of limited forward-Euler outputs and previously accepted data.

For \(M=0\), the limited forward-Euler step is TVDM under the appropriate CFL condition. Since \(\operatorname{TVM}\) is a seminorm and 12 is a convex combination of forward-Euler steps, the TVDM inequality ?? is inherited by the SSP-RKDG method. With the TVB limiter \(M>0\), the same convexity argument gives its TVBM version.

The essential point is that the classical proof uses the SSP decomposition of the time integrator. TVD/TVB stability is obtained because every stage can be interpreted as a convex combination of limited forward-Euler updates. This is a powerful framework, but it also ties the limiting strategy to the SSP structure of the RK method. Generic RK methods that do not admit such a convex forward-Euler representation are therefore outside this classical argument.

3 New Framework for non-SSP Methods↩︎

We now introduce a framework distinct from the SSP argument in Section 2, which treats a fully discrete scheme as a sequence of predictors at intermediate stages followed by a corrector at the final stage. The TVD/TVB structure is then imposed solely by modifying the numerical fluxes in the final corrector. This approach is independent of the predictors, and hence it applies universally to the standard RKDG, cRKDG, and ADER-DG methods. In the standard RKDG case, the predictors are computed via RK approximations using the usual DG operator \(\mathcal{F}\); in the cRKDG case, the predictors are computed via RK approximations using a compact cell-local operator \(\mathcal{G}\); in the ADER-DG case, the predictors are computed via a local space-time Galerkin method.

To maintain focus, we concentrate on RK formulations, discussing RKDG and cRKDG schemes in detail, while only briefly commenting on the application to the ADER-DG method in Section 3.5.

3.1 Generic predictor-corrector formulations↩︎

Consider an explicit \(s\)-stage RK method with the Butcher tableau \[\label{eq:Butcher-RK} \begin{array}{c|c} c & A \\ \hline & b \end{array} \quad c = (c_1,\cdots,c_s)^\top, \; b = (b_1,\cdots,b_s), \; A = (a_{\ell j})_{s\times s}, \; a_{\ell j}=0 \text{ if } j \ge \ell.\tag{13}\]

The standard RKDG method applies this RK integrator to the semi-discrete DG system 2 . Starting from an accepted solution \(u_h^n\), the predictor stages are \[u_h^{(\ell)} = u_h^n +\Delta t\sum_{j=1}^{\ell-1}a_{\ell j}\mathcal{F}(u_h^{(j)}), \qquad \ell=1,\ldots,s, \label{eq:rkdg-predictor}\tag{14}\] and the corresponding raw RKDG corrector is \[u_h^{n+1} = u_h^n +\Delta t\sum_{\ell=1}^s b_\ell\mathcal{F}(u_h^{(\ell)}). \label{eq:rkdg-corrector-raw}\tag{15}\] If the RK method is not SSP, then the argument in Section 2.3 does not provide a TVDM/TVBM mechanism for 14 –15 .

For the cRKDG method, we define the cell-local operator \(\mathcal{G}:\mathbb{V}_h^k\to \mathbb{V}_h^k\) by \[\int_{I_i}\mathcal{G}(u_h)v\,\mathrm{d}x = \int_{I_i}f(u_h)v_x\,\mathrm{d}x - \left[ f(u_{h,i+1/2}^{-})v_{i+1/2}^{-} - f(u_{h,i-1/2}^{+})v_{i-1/2}^{+} \right] \label{eq:G-def}\tag{16}\] for all \(v\in \mathbb{V}_h^k\). Compared with \(\mathcal{F}\), the operator \(\mathcal{G}\) uses only the inner traces from the current cell and therefore acts locally on each cell. The cRKDG method replaces \(\mathcal{F}\) by \(\mathcal{G}\) in the predictor stages 14 and keeps \(\mathcal{F}\) in the corrector 15 .

It is helpful to write both RKDG and cRKDG predictors in the unified form \[\label{eq:generic-predictor} u_h^{(\ell)} = u_h^n +\Delta t\sum_{j=1}^{\ell-1}a_{\ell j}\mathcal{H}(u_h^{(j)}), \qquad \ell=1,\ldots,s, \qquad \mathcal{H}\in\{\mathcal{F},\mathcal{G}\}.\tag{17}\] The trace-limited corrector to be introduced is the same for both choices.

3.2 Trace limiting and modified corrector operator↩︎

The polynomial limiter \(\Lambda\) defined in Section 2.2 returns a DG polynomial and preserves cell averages. The new limiter \(\Lambda_{\mathrm{tr}}\) introduced here has a different purpose. It modifies only the traces used in the final numerical flux and does not reconstruct a DG polynomial.

We denote by \(\mathbb{T}_h^k = \{ v|_{\partial \mathcal{T}_h} : v \in \mathbb{V}_h^k \}\) the space of double-valued traces for functions in \(\mathbb{V}_h^k\), where \(\partial \mathcal{T}_h = \bigcup_{i=1}^N \{ \partial I_i \}\) represents the collection of all element boundaries. Let \(v \in \mathbb{V}_h^k\) be a DG function providing the raw traces to be modified, and let \(\bar{q} \in \mathbb{V}_h^0\) be a piecewise constant function prescribing the reference cell averages. The trace limiter is defined as the mapping \(\Lambda_{\mathrm{tr}}: \mathbb{V}_h^k \times \mathbb{V}_h^0 \to \mathbb{T}_h^k\), where the output \(\widetilde{v} = \Lambda_{\mathrm{tr}}(v, \bar{q}) = \{ \widetilde{v}_{i-1/2}^+, \widetilde{v}_{i+1/2}^- \}_{i=1}^N\) denotes the limited traces on \(\partial \mathcal{T}_h\) such that \[\tag{18} \begin{align} \widetilde{v}_{i+1/2}^{-} &= \overline{q}_i + \widetilde{m} \left( v_{i+1/2}^{-}-\overline{q}_i, \Delta_+\overline{q}_i, \Delta_-\overline{q}_i; M,h_i \right), \tag{19} \\ \widetilde{v}_{i-1/2}^{+} &= \overline{q}_i - \widetilde{m} \left( \overline{q}_i-v_{i-1/2}^{+}, \Delta_+\overline{q}_i, \Delta_-\overline{q}_i; M,h_i \right). \tag{20} \end{align}\] Note that 18 generalizes 8 by allowing \(\overline{q}_i \neq \overline{v}_{i}\). The raw traces and reference cell averages can be completely unrelated.

We now generalize the DG operator in 3 by separating the volume and facial contributions, which yields \(\mathcal{F}(\cdot;\cdot):\mathbb{V}_h^k\times \mathbb{T}_h^k\to \mathbb{V}_h^k\). Given \(v\in \mathbb{V}_h^k\) and limited traces \(\widetilde{v} \in \mathbb{T}_h^k\), we define \(\mathcal{F}(v;\widetilde{v})\in \mathbb{V}_h^k\) by setting \[\begin{align} \int_{I_i}\mathcal{F}(v;\widetilde{v})\phi\,\mathrm{d}x =\int_{I_i}f(v)\phi_x\,\mathrm{d}x - \left[ \widehat f(\widetilde{v}_{i+\frac{1}{2}}^{-},\widetilde{v}_{i+\frac{1}{2}}^{+}) \phi_{i+\frac{1}{2}}^{-} - \widehat f(\widetilde{v}_{i-\frac{1}{2}}^{-},\widetilde{v}_{i-\frac{1}{2}}^{+}) \phi_{i-\frac{1}{2}}^{+} \right] \end{align} \label{eq:F-two-argument}\tag{21}\] for all \(\phi\in \mathbb{V}_h^k\). The volume term is evaluated using the input raw polynomial \(v\), while only the numerical fluxes are evaluated using the limited traces \(\widetilde{v}\). If the limited traces coincide with the original traces of \(v\), then \(\mathcal{F}(v;\widetilde{v})=\mathcal{F}(v)\).

3.3 Proposed trace-limited scheme↩︎

We now define the proposed method. Let \(u_h^n\in \mathbb{V}_h^k\) be the accepted DG solution at time \(t^n\), and let \(\overline{u}_h^n \in \mathbb{V}_h^0\) be the piecewise constant function defined by the cell averages. First, we compute the predictor stages via 17 , with either the RKDG or cRKDG predictors. Next, apply the trace limiter to each predictor stage using the initial cell averages \(\overline{u}_h^n\) as the reference: \[\widetilde{u}_h^{(\ell)} = \Lambda_{\mathrm{tr}} \left( u_h^{(\ell)}; \overline{u}_h^n \right), \qquad \ell=1,\ldots,s. \label{eq:proposed-trace-limiting}\tag{22}\] Then, the corrected candidate at the new time level is defined by \[w_h^{n+1} = u_h^n +\Delta t\sum_{\ell=1}^s b_\ell\mathcal{F}(u_h^{(\ell)};\widetilde{u}_h^{(\ell)}). \label{eq:proposed-corrector}\tag{23}\] Finally, the accepted solution is obtained by the polynomial limiter \[u_h^{n+1} = \Lambda(w_h^{n+1}). \label{eq:proposed-post-limiter}\tag{24}\] Equations 17 , 22 –24 form the proposed trace-limited RKDG/cRKDG scheme.

The proposed scheme separates the predictors from the limiting mechanism. In 17 , the choice \(\mathcal{H}=\mathcal{F}\) gives a trace-limited RKDG method, while \(\mathcal{H}=\mathcal{G}\) gives a trace-limited cRKDG method. The final corrector 23 and the final polynomial limiter 24 are identical in both cases.

Throughout the analysis below, we assume that the corrector weights satisfy \[b_\ell\ge 0, \qquad \sum_{\ell=1}^s b_\ell=1. \label{eq:positive-b-weights}\tag{25}\] This condition is much weaker than requiring the full RK method to be SSP. It is satisfied by many commonly used RK integrators, such as Heun’s third-order method and the classical fourth-order RK method discussed in the following subsection. However, some specially constructed RK time integrators may not satisfy this condition, such as the counterexample in [7] for demonstrating the necessity of SSP time integrators.

3.4 Concrete cRKDG schemes↩︎

We now list concrete cRKDG schemes used in the numerical tests. The RKDG versions are obtained by replacing \(\mathcal{G}\) by \(\mathcal{F}\).

Heun’s third-order method: \[\begin{align} u_h^{(1)}&=u^{n}_{h},\;\; u_h^{(2)} = u_h^n+\frac{\Delta t}{3}\mathcal{G}(u_h^{(1)}),\;\; u_h^{(3)} = u_h^n+\frac{2\Delta t}{3}\mathcal{G}(u_h^{(2)}), \\ \widetilde{u}_h^{(\ell)} &= \Lambda_{\mathrm{tr}} \left(u_h^{(\ell)};\overline{u}_h^n\right),\quad \ell=1,3, \\ w_h^{n+1} &= u_h^n +\Delta t \left[ \frac{1}{4}\mathcal{F}(u_h^{(1)};\widetilde{u}_h^{(1)}) + \frac{3}{4}\mathcal{F}(u_h^{(3)};\widetilde{u}_h^{(3)}) \right], \\ u_h^{n+1} &= \Lambda(w_h^{n+1}). \end{align} \label{eq:crkdg3-limited}\tag{26}\]

Classic fourth-order method: \[\begin{align} u_h^{(1)}&=u_{h}^{n},\; u_h^{(2)} = u_h^n+\frac{\Delta t}{2}\mathcal{G}(u_h^{(1)}), \; u_h^{(3)} = u_h^n+\frac{\Delta t}{2}\mathcal{G}(u_h^{(2)}), \; u_h^{(4)} = u_h^n+\Delta t\,\mathcal{G}(u_h^{(3)}), \\ \widetilde{u}_h^{(j)} &= \Lambda_{\mathrm{tr}} \left(u_h^{(j)};\overline{u}_h^n\right), \qquad j=1,2,3,4, \\ w_h^{n+1} &= u_h^n +\Delta t \left[ \frac{1}{6}\mathcal{F}(u_h^{(1)};\widetilde{u}_h^{(1)}) + \frac{1}{3}\mathcal{F}(u_h^{(2)};\widetilde{u}_h^{(2)}) + \frac{1}{3}\mathcal{F}(u_h^{(3)};\widetilde{u}_h^{(3)}) + \frac{1}{6}\mathcal{F}(u_h^{(4)};\widetilde{u}_h^{(4)}) \right], \\ u_h^{n+1} &= \Lambda(w_h^{n+1}). \end{align} \label{eq:crkdg4-limited}\tag{27}\]

3.5 Applications to ADER-DG methods↩︎

Although this paper primarily focuses on RKDG and cRKDG methods, the proposed predictor-corrector framework naturally accommodates the ADER-DG method, which is derived from a space-time weak formulation of the conservation law. Below, we discuss one version of this method that utilizes a local predictor [16], [18], [34].

The first step is to compute a predictor using a local space-time Galerkin method for 1 : find \(z_h \in \mathbb{P}^{k'}(I_i \times [t^n,t^{n+1}])\), such that for all \(w_h \in \mathbb{P}^{k'}(I_i\times[t^n,t^{n+1}])\), \[\label{eq:ADER-I} \int_{I_i} z_h(x,t^{n+1}) w_h(x,t^{n+1}) - u_h^n w_h(x,t^{n}) \, \mathrm{d} x +\int_{t^n}^{t^{n+1}}\int_{I_i} -z_h(w_h)_t + f(z_h)_x w_h \, \mathrm{d} x \, \mathrm{d} t = 0.\tag{28}\] Then, the corrector \(u_h^{n+1} \in \mathbb{V}_h^k\) is updated via \[\label{eq:ADER-II} u_h^{n+1} = u_h^n + \int_{t^n}^{t^{n+1}} \mathcal{F}(z_h(x,t)) \, \mathrm{d} t.\tag{29}\]

Applying a quadrature rule on \([t^n,t^{n+1}]\) with nodes \(\{t^{n,\ell}\}\) and weights \(\{b_\ell\}\), one can rewrite 29 in a form similar to 15 , where \(u_h^{(\ell)}(x) = z_h(x,t^{n,\ell})\). Consequently, 28 , the relation \(u_h^{(\ell)}(x) = z_h(x,t^{n,\ell})\), and 22 –24 constitute the trace-limited ADER-DG scheme. This scheme shares similar theoretical properties with its RKDG and cRKDG counterparts, and the details are omitted.

4 Properties of Trace-Limited Schemes↩︎

In this section, we establish the main properties of the trace-limited scheme 17 , 22 –24 for the 1D scalar equation 1 with periodic or compactly supported boundaries. Throughout the section, we assume a uniform mesh with cell size \(h\), and write \(\lambda={\Delta t}/{h}.\)

4.1 Conservation↩︎

We first note that the trace-limited corrector preserves the conservative form of the DG method. Although the trace limiter modifies the left and right traces entering the numerical flux, the numerical flux at each interface is still single-valued. Therefore the flux leaving one cell is the same as the flux entering the neighboring cell. This observation leads to the following proposition on conservation, whose detailed proof is omitted.

4.2 TVDM and TVBM properties↩︎

The key of the proposed method is that the TVD mechanism depends only on the final cell-average update, not on how the predictor stages are computed. The following lemma makes this precise. It shows that a forward-Euler-type update is TVDM even if the raw traces entering the limiter are arbitrary, provided that the limiter uses the old cell averages as reference values.

We now apply Lemma 2 to the full trace-limited RK corrector. The result uses only the final corrector 23 and it does not depend on whether the predictor stages were computed with \(\mathcal{F}\) or \(\mathcal{G}\).

4.3 Non-activation in smooth regions↩︎

Finally, we address the accuracy of the limiting procedure. Under a mild assumption on the predictor stages and a CFL condition, we show that, for sufficiently small \(h\), the trace limiter remains inactive in smooth monotone regions away from critical points, and also remains inactive near smooth extrema when the TVB constant \(M\) is sufficiently large. Consequently, the limiting procedure preserves the formal order of accuracy of the underlying RKDG or cRKDG schemes.

Recall that \(c_\ell\) denotes the Butcher node for stage \(\ell\), with \(c_\ell=\sum_{j=1}^{\ell-1} a_{\ell j}\).

Theorem 2 should be read as an accuracy-preservation statement for the limiter, not as a standalone convergence theorem for the full RKDG or cRKDG method. It says that if the predictor stages already have the expected local accuracy in smooth regions, then the trace limiter is inactive there for sufficiently fine meshes. The post-step polynomial limiter \(\Lambda\) satisfies the analogous standard TVB non-activation property when applied to the smooth candidate \(w_h^{n+1}\) [2].

5 Extensions to Systems and Multi-Dimensions↩︎

While previous sections focused on 1D scalar conservation laws, we now extend the method to systems and multidimensional problems. The time discretization and trace-limited corrector structure remain unchanged; what changes is only the application of the trace limiter \(\Lambda_{\mathrm{tr}}\) to vector-valued or multidimensional DG polynomials. While our focus here is on the 2D case, this extension applies naturally to higher dimensions.

5.1 1D systems↩︎

Consider the 1D system \[\mathbf{u}_t+\mathbf{f}(\mathbf{u})_x=0, \qquad \mathbf{u}\in\mathbb{R}^m. \label{eq:system-1d}\tag{39}\] The trace-limited RKDG/cRKDG scheme has the same form as that in Section 3.3, with all scalar quantities replaced by vectors: \[\label{eq:system-scheme-1d} \begin{align} \mathbf{u}_h^{(\ell)} =& \mathbf{u}_h^n +\Delta t\sum_{j=1}^{\ell-1}a_{\ell j}\mathcal{H}(\mathbf{u}_h^{(j)}),\qquad \ell=1,\ldots,s, \qquad \mathcal{H}\in\{\mathcal{F},\mathcal{G}\},\\ \widetilde{\mathbf{u}}_h^{(\ell)} =& \Lambda_{\mathrm{tr}} \left(\mathbf{u}_h^{(\ell)};\overline{\mathbf{u}}_h^n\right), \qquad \ell=1,\ldots,s.\\ \mathbf{w}_h^{n+1} =& \mathbf{u}_h^n +\Delta t\sum_{\ell=1}^s b_\ell \mathcal{F}(\mathbf{u}_h^{(\ell)};\widetilde{\mathbf{u}}_h^{(\ell)}),\\ \mathbf{u}_h^{n+1}=&\Lambda(\mathbf{w}_h^{n+1}). \end{align}\tag{40}\]

We now define the trace limiter \(\Lambda_{\mathrm{tr}}:[\mathbb{V}_h^k]^m\times[\mathbb{V}_h^0]^m \to [\mathbb{T}_h^k]^m\). Throughout this section, for a given pair \((\mathbf{v},\overline{\mathbf{q}})\in [\mathbb{V}_h^k]^m\times[\mathbb{V}_h^0]^m\), we write \[\widetilde{\mathbf{v}} = \Lambda_{\mathrm{tr}}(\mathbf{v};\overline{\mathbf{q}})\] for the limited trace function. On a single mesh cell \(I_i\), the quantities \(\mathbf{v}_{i+1/2}^-\) and \(\mathbf{v}_{i-1/2}^+\) denote the raw traces to be limited, while \(\overline{\mathbf{q}}_i\) denotes the reference cell average used by the limiter. Define the trace deviations from the reference average by \[\boldsymbol{\delta}_{i,+} = \mathbf{v}_{i+1/2}^{-}-\overline{\mathbf{q}}_i, \qquad \boldsymbol{\delta}_{i,-} = \overline{\mathbf{q}}_i-\mathbf{v}_{i-1/2}^{+}. \label{eq:system-raw-deviations-1d}\tag{41}\] The neighboring mean differences are defined by \[\Delta_+ \overline{\mathbf{q}}_{i}= \overline{\mathbf{q}}_{i+1}-\overline{\mathbf{q}}_{i}, \qquad \Delta_- \overline{\mathbf{q}}_{i}= \overline{\mathbf{q}}_{i}-\overline{\mathbf{q}}_{i-1}.\]

For systems, the limiting is performed in local characteristic variables. We need to use the eigen-decomposition \(A_i =\mathbf{f}'(\overline{\mathbf{q}}_i) = R_i D_i R_i^{-1}\), where \(D_i\) is a diagonal matrix of eigenvalues and the columns of \(R_i\) are the corresponding right eigenvectors. These quantities are transformed into characteristic variables: \[\begin{align} \boldsymbol{\alpha}_{i,+} =R_i^{-1}\boldsymbol{\delta}_{i,+}, \quad \boldsymbol{\alpha}_{i,-} =R_i^{-1}\boldsymbol{\delta}_{i,-}, \quad \boldsymbol{\alpha}_{i}^{+} =R_i^{-1}\Delta_+ \overline{\mathbf{q}}_{i}, \quad \boldsymbol{\alpha}_{i}^{-} =R_i^{-1}\Delta_-\overline{\mathbf{q}}_{i}. \end{align} \label{eq:system-char-differences-1d}\tag{42}\] The scalar TVB-modified minmod function is then applied componentwise: \[\begin{align} \widetilde{\boldsymbol{\alpha}}_{i,+} &= \widetilde{m} \left( \boldsymbol{\alpha}_{i,+}, \boldsymbol{\alpha}_{i}^{+}, \boldsymbol{\alpha}_{i}^{-};\mathbf{M},h_i \right),\qquad \widetilde{\boldsymbol{\alpha}}_{i,-} = \widetilde{m} \left( \boldsymbol{\alpha}_{i,-}, \boldsymbol{\alpha}_{i}^{+}, \boldsymbol{\alpha}_{i}^{-};\mathbf{M},h_i \right). \end{align} \label{eq:system-char-minmod-1d}\tag{43}\] The TVB parameter \(\mathbf{M}\in\mathbb{R}^{m}\) may also be chosen componentwise in characteristic variables. The limited traces are defined by transforming back the modified values to conservative variables and adding them to the cell averages: \[\widetilde{\mathbf{v}}_{i+1/2}^{-} = \overline{\mathbf{q}}_i+R_i\widetilde{\boldsymbol{\alpha}}_{i,+}, \qquad \widetilde{\mathbf{v}}_{i-1/2}^{+} = \overline{\mathbf{q}}_i-R_i\widetilde{\boldsymbol{\alpha}}_{i,-}. \label{eq:system-limited-traces-1d}\tag{44}\] This completes the definition of trace limiting \(\widetilde{\mathbf{v}} = \Lambda_\mathrm{tr}(\mathbf{v}, \overline{\mathbf{q}})\).

The polynomial limiter \(\Lambda\) for systems is defined in the standard way. When limiting a polynomial \(\mathbf{v}\), the reference average in the trace limiter is set as its own cell average, i.e. \(\overline{\mathbf{q}}=\overline{\mathbf{v}}\). A cell \(I_i\) is called troubled if at least one of its traces is modified by \(\Lambda_{\mathrm{tr}}(\mathbf{v};\overline{\mathbf{v}})\). On nontroubled cells, \(\Lambda\) leaves the polynomial \(\mathbf{v}|_{I_i}\) unchanged. On a troubled cell \(I_i\), \(\Lambda\) replaces \(\mathbf{v}|_{I_i}\) with a polynomial in \([\mathbb{P}^{\min(k,2)}(I_i)]^m\) that preserves the cell average \(\overline{\mathbf{v}}_i\) and matches the limited interface traces \(\widetilde{\mathbf{v}}_{i+1/2}^-\) and \(\widetilde{\mathbf{v}}_{i-1/2}^+\).

5.2 2D scalar equations↩︎

We next consider a scalar conservation law in 2D, \[u_t+f(u)_x+g(u)_y=0,\qquad (x,y)\in\Omega,\quad t>0 . \label{eq:scalar-2d}\tag{45}\] Let \(\mathcal{T}_h\) be a Cartesian mesh with cells \(K_{ij}=I_i\times J_j=[x_{i-1/2},x_{i+1/2}] \times[y_{j-1/2},y_{j+1/2}]\). We set \(h_i^x=x_{i+1/2}-x_{i-1/2}\), \(h_j^y=y_{j+1/2}-y_{j-1/2}\), and introduce the reference coordinates \(\xi=2(x-x_i)/h_i^x\), \(\eta=2(y-y_j)/h_j^y\), with \((\xi,\eta)\in[-1,1]^2\), on each cell \(K_{ij}\). The DG space is \(\mathbb{V}_h^k=\{v\in L^2(\Omega):v|_{K_{ij}}\in \mathbb{Z}^k(K_{ij}) \;\text{for all } i,j\}\), where \(\mathbb{Z}^k(K_{ij})\) denotes either the total-degree space \(\mathbb{P}^k(K_{ij})\) or the tensor-product space \(\mathbb{Q}^k(K_{ij})\).

We now define the trace limiter \(\Lambda_{\mathrm{tr}}:\mathbb{V}_h^k\times\mathbb{V}_h^0\to\mathbb{T}_h^k\), where \(\mathbb{T}_h^k=\{v|_{\partial\mathcal{T}_h}:v\in\mathbb{V}_h^k\}\) and \(\partial\mathcal{T}_h=\{\partial K:K\in\mathcal{T}_h\}\). For a given pair \((v,\overline{q})\in \mathbb{V}_h^k\times \mathbb{V}_h^0\), we write \(\widetilde{v}=\Lambda_{\mathrm{tr}}(v;\overline{q})\) for the limited trace function. Let the affine part (the \(\mathbb{P}^1\)-component) of \(v\) on \(K_{ij}\) be \[\Pi_1 v|_{K_{ij}} = \overline{v}_{ij} + v^x_{ij}\xi + v^y_{ij}\eta . \label{eq:P1-projection-2d}\tag{46}\] The edge-average deviations of this affine part relative to \(\overline{q}_{ij}\) are \[\delta_R=(\overline{v}_{ij}+v^x_{ij})-\overline{q}_{ij},\; \delta_L=\overline{q}_{ij}-(\overline{v}_{ij}-v^x_{ij}),\; \delta_T=(\overline{v}_{ij}+v^y_{ij})-\overline{q}_{ij},\; \delta_B=\overline{q}_{ij}-(\overline{v}_{ij}-v^y_{ij}),\] corresponding to the right, left, top, and bottom edges \(\partial K_{ij}^{E}, E\in\{R,L,T,B\}\), respectively. The limited deviations are then defined by \[\tag{47} \begin{align} \widetilde{\delta}_E &= \widetilde{m} \left( \delta_E, \overline{q}_{i+1,j}-\overline{q}_{ij}, \overline{q}_{ij}-\overline{q}_{i-1,j}; M,h_i^x \right), \qquad E\in\{R,L\}, \tag{48} \\ \widetilde{\delta}_E &= \widetilde{m} \left( \delta_E, \overline{q}_{i,j+1}-\overline{q}_{ij}, \overline{q}_{ij}-\overline{q}_{i,j-1}; M,h_j^y \right), \qquad E\in\{T,B\}. \tag{49} \end{align}\]

If none of the four deviations is modified, then \(\Lambda_{\mathrm{tr}}(v;\overline{q})\) returns the original inner traces of \(v\) on \(\partial K_{ij}\). Otherwise, construct the limited affine function \[\widetilde{v}_{ij}(\xi,\eta) = \overline{q}_{ij} + \frac{1}{2}(\widetilde{\delta}_R+\widetilde{\delta}_L)\xi + \frac{1}{2}(\widetilde{\delta}_T+\widetilde{\delta}_B)\eta . \label{eq:modaffine}\tag{50}\] For cells with active limiting, the trace limiter returns four inner edge traces of 50 : \[\Lambda_{\mathrm{tr}}(v;\overline{q})|_{\partial K_{ij}^{E}} = \widetilde{v}_{ij}|_{\partial K_{ij}^{E}}, \qquad E\in\{R,L,T,B\}. \label{eq:Lambda-trace-2d-scalar}\tag{51}\]

The polynomial limiter \(\Lambda:\mathbb{V}_h^k\to\mathbb{V}_h^k\) is the classical Cockburn–Shu limiter in multiple dimensions on Cartesian meshes [5]. In the above construction, it corresponds to taking \(\overline{q}=\overline{v}\). Thus, on a cell \(K_{ij}\), \((\Lambda v)|_{K_{ij}}=v|_{K_{ij}}\) if the limiting in 47 is inactive, while \((\Lambda v)|_{K_{ij}}=\widetilde{v}_{ij}\) otherwise, where \(\widetilde{v}_{ij}\) is defined by 50 with \(\overline{q}_{ij}=\overline{v}_{ij}\).

The next proposition shows that the TVB-modified minmod limiter is idempotent on active cells; re-applying it to the resulting edge averages leaves them unchanged.

5.3 2D systems↩︎

For a 2D system \[\mathbf{u}_t+\mathbf{f}(\mathbf{u})_x+\mathbf{g}(\mathbf{u})_y=0, \qquad \mathbf{u}\in\mathbb{R}^m, \label{eq:system-2d}\tag{55}\] the limiters are applied in local characteristic variables, separately in the two coordinate directions. The construction combines the characteristic limiting in Section 5.1 with the edge-average affine reconstruction in Section 5.2.

We now define the trace limiter \(\Lambda_{\mathrm{tr}}:[\mathbb{V}_h^k]^m\times[\mathbb{V}_h^0]^m \to[\mathbb{T}_h^k]^m.\) For \(\mathbf{v}\in[\mathbb{V}_h^k]^m\) and \(\overline{\mathbf{q}}\in[\mathbb{V}_h^0]^m\), write \(\widetilde{\mathbf{v}}=\Lambda_{\mathrm{tr}}(\mathbf{v};\overline{\mathbf{q}})\). On a cell \(K_{ij}\), consider the eigen-decompositions \[A^x_{ij} = \mathbf{f}'(\overline{\mathbf{q}}_{ij}) = R^x_{ij}D^x_{ij}(R^x_{ij})^{-1}, \qquad A^y_{ij} = \mathbf{g}'(\overline{\mathbf{q}}_{ij}) = R^y_{ij}D^y_{ij}(R^y_{ij})^{-1},\] where the columns of \(R^x_{ij}\) and \(R^y_{ij}\) are the corresponding right eigenvectors. Write the affine part of \(\mathbf{v}\) on \(K_{ij}\) as \[\Pi_1\mathbf{v}|_{K_{ij}} = \overline{\mathbf{v}}_{ij} + \mathbf{v}^x_{ij}\xi + \mathbf{v}^y_{ij}\eta . \label{eq:P1-projection-2d-system}\tag{56}\] The edge-average deviations relative to \(\overline{\mathbf{q}}_{ij}\) are \[\begin{align} \boldsymbol{\delta}_R&=(\overline{\mathbf{v}}_{ij}+\mathbf{v}^x_{ij}) -\overline{\mathbf{q}}_{ij}, &\quad \boldsymbol{\delta}_L&=\overline{\mathbf{q}}_{ij} -(\overline{\mathbf{v}}_{ij}-\mathbf{v}^x_{ij}),\\ \boldsymbol{\delta}_T&=(\overline{\mathbf{v}}_{ij}+\mathbf{v}^y_{ij}) -\overline{\mathbf{q}}_{ij}, &\quad \boldsymbol{\delta}_B&=\overline{\mathbf{q}}_{ij} -(\overline{\mathbf{v}}_{ij}-\mathbf{v}^y_{ij}). \end{align}\] The neighboring mean differences in \(x\)- and \(y\)-directions are defined by \[\Delta_\pm^x\overline{\mathbf{q}}_{ij} = \pm(\overline{\mathbf{q}}_{i\pm1,j}-\overline{\mathbf{q}}_{ij}), \qquad \Delta_\pm^y\overline{\mathbf{q}}_{ij} = \pm(\overline{\mathbf{q}}_{i,j\pm1}-\overline{\mathbf{q}}_{ij}).\]

For \(d=x,y\), set \(\mathcal{E}_x=\{R,L\}\) and \(\mathcal{E}_y=\{T,B\}\). For each edge \(E\in\mathcal{E}_d\), project the corresponding edge deviation and neighboring mean differences into the characteristic variables of \(A^d_{ij}\): \[\boldsymbol{\alpha}_E^d=(R^d_{ij})^{-1}\boldsymbol{\delta}_E, \qquad \boldsymbol{\alpha}_\pm^d=(R^d_{ij})^{-1} \Delta_\pm^d\overline{\mathbf{q}}_{ij}.\] The scalar TVB-modified minmod function is then applied componentwise, and the result is transformed back to conservative variables: \[\widetilde{\boldsymbol{\delta}}_E = R^d_{ij}\, \widetilde{m} \left( \boldsymbol{\alpha}_E^d, \boldsymbol{\alpha}_+^d, \boldsymbol{\alpha}_-^d; \mathbf{M},h^d_{i/j} \right), \quad E\in\mathcal{E}_d,\quad d=x,y, \quad h_{i/j}^d = h_i^x\text{ or } h_j^y.\]

If the componentwise characteristic limiting does not modify any of the four edge deviations, then \(\Lambda_{\mathrm{tr}}(\mathbf{v};\overline{\mathbf{q}})\) returns the original inner traces of \(\mathbf{v}\). Otherwise, construct \[\widetilde{\mathbf{v}}_{ij}(\xi,\eta) = \overline{\mathbf{q}}_{ij} + \frac{1}{2}(\widetilde{\boldsymbol{\delta}}_R+ \widetilde{\boldsymbol{\delta}}_L)\xi + \frac{1}{2}(\widetilde{\boldsymbol{\delta}}_T+ \widetilde{\boldsymbol{\delta}}_B)\eta. \label{eq:system-2d-affine-reconstruction}\tag{57}\] On an active cell, the trace limiter returns the inner edge traces of \(\widetilde{\mathbf{v}}_{ij}\).

The polynomial limiter \(\Lambda:[\mathbb{V}_h^k]^m\to[\mathbb{V}_h^k]^m\) leaves \(\mathbf{v}|_{K_{ij}}\) unchanged on nontroubled cells, while it replaces \(\mathbf{v}|_{K_{ij}}\) by the limited affine reconstruction \(\widetilde{\mathbf{v}}_{ij}\) with \(\overline{\mathbf{q}}=\overline{\mathbf{v}}\) on troubled cells. This retrieves the usual characteristic-wise Cockburn–Shu limiter for systems on Cartesian meshes [5].

6 Numerical Tests↩︎

We test the proposed trace-limited cRKDG schemes on a suite of scalar conservation laws and the Euler equations, and compare the results with those obtained using classical SSP-RKDG schemes of the same order. To save space, we present only the results for third- and fourth-order schemes in the accuracy tests, and for third-order schemes in the shock-capturing tests. Besides 12 , 26 and 27 , we use the Spiteri–Ruuth scheme in [36] for the fourth-order SSP-RKDG method. All schemes use the local Lax–Friedrichs numerical flux, with the dissipation parameter set to \(\alpha=1.05\) times the local wave speed. We set \(M=50\) uniformly in all TVB tests, following the default choice in [5]. The CFL numbers are set to \(0.05\) for the accuracy tests and \(0.15\) for all shock-capturing tests. The orders of spatial accuracy are chosen to be consistent with the corresponding temporal orders.

6.1 Scalar conservation laws↩︎

6.1.1 Linear advection↩︎

We consider the linear advection equation \(u_t + u_x = 0\) on the periodic domain \([0,2\pi]\).

6.1.1.1 Accuracy test

We take the smooth initial data \(u(x,0) = \sin(x)\), so that the exact solution is \(u(x,t) = \sin(x-t)\). The final time is \(T=1\). Table 1 reports the \(L^2\) errors and convergence orders for the schemes under the TVD (\(M=0\)) and TVB (\(M=50\)) settings, respectively. From the table, we observe optimal accuracy with TVB limiting and order degeneration in the TVD setting for all schemes.

6.1.1.2 Shock-capturing test

We solve the same equation with the discontinuous square-wave initial condition \(u(x,0) = 1\) for \(\pi/2 \le x \le 3\pi/2\), and \(u(x,0) = 0\) otherwise, using \(N=80\) cells. The final time is \(T=2\pi\), so the exact solution travels one period. Figures 1 (a) and 1 (b) compare the cell averages of \(\mathbb{P}^2\) cRK3 and \(\mathbb{P}^2\) SSP-RK3 solutions against the exact solution under the TVD and TVB limiters, respectively. From the figures we can observe almost identical TVDM behaviors of both schemes, and comparable results in the TVB setting.

6.1.2 Burgers equation↩︎

We consider the Burgers equation \(u_t + ({u^2}/{2})_x = 0\) on the periodic domain \([0,2\pi]\).

6.1.2.1 Accuracy test

We take the smooth initial data \(u(x,0) = {1}/{2}+\sin(x)\). The exact solution is obtained by the method of characteristics and Newton iteration. The final time \(T=0.2\) is before shock formation, so the solution remains smooth throughout the computation. Table 2 reports the \(L^2\) errors and convergence orders. As in the linear case, we observe optimal convergence orders under the TVB limiting and order degeneration under the TVD limiting.

6.1.2.2 Shock-capturing test

We solve the problem with the discontinuous square-wave initial condition \(u(x,0) = 1\) for \(\pi/2 \le x \le 3\pi/2\), and \(u(x,0) = -1\) otherwise, using \(N=80\) cells and a final time of \(T=1\). At \(x=\pi/2\) the jump from \(-1\) to \(1\) produces a rarefaction fan, while at \(x=3\pi/2\) the jump from \(1\) to \(-1\) produces a stationary shock. This test simultaneously examines the resolution of a rarefaction wave and the sharpness of a standing shock. Figures 2 (a) and 2 (b) compare the cell averages of \(\mathbb{P}^2\) cRK3 and \(\mathbb{P}^2\) SSP-RK3 solutions under the TVD and TVB limiters, respectively.

6.1.3 Buckley–Leverett equation↩︎

We consider the Buckley–Leverett equation \(u_t + f(u)_x = 0\), where the flux function is given by \(f(u) = 4u^2 / (4u^2 + (1-u)^2)\). The flux \(f\) is nonconvex in this problem. We test two Riemann problems. For both tests, the domain is \([-0.5,\,0.5]\) with \(N=80\) cells, the final time is \(T=1\), and inflow boundary conditions are imposed using the corresponding constant states.

6.1.3.1 Shock–rarefaction–shock structure

The initial data are \(u(x,0) = 2\) for \(x < 0\), and \(u(x,0) = -2\) for \(x \ge 0\). The entropy solution is a composite shock–rarefaction–shock pattern. Figures 3 (a) and 3 (b) show the results.

6.1.3.2 Two-shock structure

The initial data are \(u(x,0) = -3\) for \(x < 0\), and \(u(x,0) = 3\) for \(x \ge 0\). The entropy solution consists of two shocks separated by a constant intermediate state \(u=0\). Figures 4 (a) and 4 (b) show the results.

6.1.4 2D Burgers equation↩︎

We consider the 2D Burgers equation \(u_t + ({u^2}/{2})_x + ({u^2}/{2})_y = 0\) on the periodic domain \([0,2\pi]^2\). The initial condition is taken as \(u(x,y,0)=1/2+\sin(x+y)\). The corresponding shock formation time is \(T_s=1/2.\)

6.1.4.1 Accuracy test

The final time \(T=0.2\) is chosen before shock formation. Table 3 reports the \(L^2\) errors and convergence orders.

6.1.4.2 Shock-capturing test

The final time \(T=0.7\) is after shock formation. Figures 5 (a) and 5 (b) show the cell averages of the solutions along the horizontal line \(y=0.5\).

6.2 Euler equations↩︎

We now test the schemes on the 1D Euler equations, \[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ m \\ E \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} m \\ m^2/\rho + p \\ (E+p)\,m/\rho \end{pmatrix} = 0,\] closed by the ideal-gas equation of state \(p = (\gamma-1)\bigl(E - m^2/(2\rho)\bigr)\). Unless stated otherwise, \(\gamma = 1.4\).

6.2.1 Accuracy test↩︎

We consider \(\gamma = 3\) with the isentropic initial data \(\rho = 1 + \sin(x)/2\), \(u = {1}/{2}+\sin(x)\), and \(p = (1 + \sin(x)/2)^3\) on the periodic domain \([0,2\pi]\). The system reduces to two decoupled Burgers equations for the variables \(r_1 = u + \sqrt{3}\,\rho\) and \(r_2 = u - \sqrt{3}\,\rho\) before shock formation. An exact solution is therefore available via the method of characteristics. The final time is \(T=0.1\). Table 4 reports the \(L^2\) errors in energy (\(E\)) and the corresponding convergence orders for all schemes under both limiter settings.

6.2.2 Blast wave problem↩︎

The blast wave problem of Woodward and Colella [37] models the interaction of strong shocks in a closed tube. The initial data are \((\rho,u,p) = (1, 0, 1000)\) for \(x < 0.1\), \((\rho,u,p) = (1, 0, 0.01)\) for \(0.1 \le x < 0.9\), and \((\rho,u,p) = (1, 0, 100)\) for \(x \ge 0.9\), on the domain \([0,1]\) with reflecting boundary conditions on both ends. The test is run with \(N=400\) cells until \(T=0.038\). There is no closed-form solution for this problem. We use a WENO5 finite-difference solution on \(10{,}000\) cells as the reference. Figures 6 (a) and 6 (b) show the density profiles.

6.2.3 Shu–Osher problem↩︎

The Shu–Osher problem [38] is a standard test used to evaluate if a scheme can simultaneously handle a strong shock and small-scale smooth structures. The initial data are \((\rho,u,p) = (3.857143, 2.629369, 10.33333)\) for \(x < -4\), and \((\rho,u,p) = (1 + 0.2\sin(5x), 0, 1)\) for \(x \ge -4\), on the domain \([-5,5]\) with \(N=200\) cells. We impose transmissive boundaries and take the final time \(T=1.8\). Since there is no closed-form solution, we use a WENO5 finite-difference solution on \(10{,}000\) cells as the reference. Figures 7 (a) and 7 (b) show the density profiles.

We then test the schemes on the 2D Euler equations \[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ E \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho u \\ \rho u^2 + p \\ \rho u v \\ (E+p)\,u \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} \rho v \\ \rho u v \\ \rho v^2 + p \\ (E+p)\,v \end{pmatrix} = 0,\] with \(p = (\gamma-1)(E - \rho(u^2+v^2)/2)\) and \(\gamma = 1.4\). Only TVB (\(M=50\)) results are presented to save space.

6.2.4 Double Mach reflection↩︎

We consider the double Mach reflection problem [37]. A Mach 10 oblique shock initially makes a \(60^\circ\) angle with the \(x\)-axis, meeting the bottom wall at \(x=1/6\). The computational domain is \([0,4]\times[0,1]\). The post-shock state is \((\rho,u,v,p) = (8,\,33\sqrt{3}/8,\,-33/8,\,116.5)\) and the pre-shock state is \((\rho,u,v,p) = (1.4,\,0,\,0,\,1)\). Along the bottom boundary, the exact post-shock condition is imposed for \(x < 1/6\) and a reflecting wall for \(x \ge 1/6\). Along the top boundary, the exact motion of the oblique shock is prescribed. Inflow and outflow conditions are imposed on the left and right boundaries, respectively. The final time is \(T=0.2\). We run the schemes with TVB limiting on two meshes: a coarse mesh of \(480\times120\) cells and a fine mesh of \(1920\times480\) cells. Figure 8 shows that the two schemes produce very close results at both resolutions.

6.2.5 Forward-facing step↩︎

We consider the forward-facing step problem [37]. A uniform Mach 3 flow enters an L-shaped tunnel with a step of height \(0.2\) located at \(x=0.6\). The computational domain is \([0,3]\times[0,1]\) with the rectangular region \([0.6,3]\times[0,0.2]\) removed. The initial condition is \((\rho,u,v,p) = (1.4,\,3,\,0,\,1)\). Inflow conditions are imposed on the left boundary, outflow on the right, and reflecting walls on all remaining boundaries. The final time is \(T=4\). We run the schemes with TVB limiting on two meshes: a coarse mesh of \(240\times80\) cells and a fine mesh of \(960\times320\) cells. Again, the cRK3 and SSP-RK3 results are nearly identical; see Figure 9.

7 Conclusions↩︎

This paper presents a novel framework for preserving the TVD and TVB properties of the DG methods using non-SSP time discretizations. The new approach uses a trace-limiting strategy and limits only the final corrector stage, regardless of how intermediate predictor stages are computed. Based on Harten’s lemma, it can be shown that the trace limiting alone can preserve the TVDM/TVBM property of the method. For further robustness, a polynomial limiter is applied at the end of the final stage to produce non-oscillatory DG polynomials. The conservation, TVDM/TVBM, and accuracy-preserving properties are proved. The proposed method decouples the TVD/TVB limiting from the SSP restriction, streamlines the algorithmic flow for concurrent computing, retains the compactness of the cRKDG framework, and accommodates fully discrete schemes beyond the RK family. The algorithmic extension to multidimensional systems is also provided. Numerical experiments with the TVD/TVB cRKDG schemes demonstrate that this new approach achieves performance comparable to classical SSP-RKDG methods. Our future work includes detailed studies on extending the framework to ADER-DG, Lax–Wendroff DG, and implicit RKDG methods, as well as its compatibility with other structure-preserving techniques.

Acknowledgment↩︎

The authors thank Professor Chi-Wang Shu of Brown University for his helpful comments and suggestions, which greatly improved the paper. The authors report the usage of AI tools and assume responsibility for all content.

References↩︎

---

1. Department of Mathematics and Statistics, Binghamton University, Binghamton, NY 13902, USA. ()↩︎

2. Department of Mathematics, The University of Alabama, Tuscaloosa, AL 35487, USA. () The work of this author was partially supported by the NSF grant DMS-2208391.↩︎