Subcell-resolution methods for advection–reaction equations with stiff sources: non-splitting high-order-in-time integration and implicit-solver compatibility on coarse meshes

Wang, Shu, Yee and Sjögreen's subcell-resolution (SR) finite difference method (J. Comput. Phys. 231(1):190–214, 2012) captures the correct locations and jump sizes of under-resolved discontinuous fronts driven by stiff reaction terms on coarse meshes, but the authors' own text acknowledges limitations: the temporal accuracy is at most second order because of Strang time splitting, implicit methods cannot be used in the reaction step because the troubled values need to be modified explicitly, and the narrow internal structure of a reaction zone containing only one or a few grid points is not resolved — resolving the sharp reaction zone, the authors remark, still requires sufficiently many grid points in that zone. A 2015 follow-up by the same group resolved the extension to multi-species, multi-reaction kinetics, but remained a fractional-step method. Non-splitting high-order detonation integrators (the Du–Yang bound-preserving line, modified exponential Runge–Kutta, compact fourth-order IMEX) do not embed subcell resolution, and under-resolved detonation methodology targets wave location only; existing SR methods are fractional-step and modify troubled values explicitly, so no SR method is compatible with implicit or IMEX treatment of the stiff reaction step. The open core is the conjunction: a coarse-mesh correct-speed method for stiff advection–reaction laws and reactive Euler systems with multi-step kinetics that is simultaneously at least third order in time without operator splitting (with the subcell correction embedded consistently in the non-splitting integrator) and compatible with implicit or IMEX treatment of the stiff reaction step.

Background

Reacting gas flow is modeled by balance laws of the form $u_t + f(u)_x = s(u)$, where the flux $f$ describes advection and the source $s(u)$ describes chemical reaction. Detonation waves are the extreme case: chemical heat release occurs on a time scale far shorter than the advective transport scale, so $s$ is called stiff. Across a detonation front the solution changes abruptly inside a narrow reaction zone (in the classical Zel'dovich–von Neumann–Döring picture, a leading shock followed by such a zone). When the computational mesh is so coarse that only one or a few grid points lie inside the reaction zone, standard shock-capturing schemes — for example those built on WENO (weighted essentially non-oscillatory) reconstruction, i.e., high-order polynomial reconstructions weighted to avoid spurious oscillations near discontinuities — propagate the front at grossly wrong speeds, because cell averages mix burnt and unburnt states and the stiff source then acts on these unphysical mixtures. Wang, Shu, Yee and Sjögreen (2012) introduced a subcell resolution (SR) remedy for this setting: cells containing the discontinuity are flagged as troubled cells (via a WENO-based smoothness indicator), and within a troubled cell a piecewise reconstruction that respects the two distinct states on either side replaces the point values before the reaction step is applied, together with per-variable extrapolation. Combined with Strang operator splitting — alternating advection and reaction sub-solves in a half-step/full-step/half-step pattern, which is second-order accurate in time — the method reproduces correct shock and front locations and jump sizes for under-resolved detonation-type fronts on very coarse meshes. The authors acknowledge limitations of the method in their condensed CTR 2011 version: (1) only single-step, two-chemical-state kinetics were treated; (2) the temporal accuracy is at most second order due to the time splitting; (3) the narrow reaction zone is not resolved — typically one or a few grid points lie inside it, and resolving the sharp reaction zone still requires sufficiently many grid points in that zone; and (4) implicit methods cannot be used in the reaction step because the troubled values need to be modified explicitly, which makes the approach incompatible with implicit stiff ODE solvers, otherwise the natural tool for stiff kinetics. The 2015 follow-up by the same group (Wang, Shu, Yee, Kotov, Sjögreen) removed limitation (1) by extending SR to multispecies, multireaction inviscid reactive flows in one and two dimensions, but remained a fractional-step method with the troubled values modified explicitly in the reaction step. Related but distinct lines that address pieces of the remaining difficulties without SR include: bound-preserving, sign-preserving non-splitting DG and finite difference methods of Du, Wang, Qian, Yang and Du–Yang (DG denotes the discontinuous Galerkin method, where the solution is represented by piecewise polynomials with weakly enforced continuity; bound-preserving means discrete mass fractions provably stay in $[0,1]$); modified exponential Runge–Kutta DG methods of Huang and Shu (Runge–Kutta, RK, being the standard family of one-step multi-stage explicit or implicit time integrators); random-projection and split-random-time-stepping correct-speed methods; and IMEX (implicit–explicit) schemes, which treat the advection part explicitly and the stiff reaction part implicitly. None of these embeds the SR troubled-point correction, and the field's under-resolved detonation methodology explicitly frames its goal as predicting wave locations without resolving internal zone structure.

Problem Statement

Construct a subcell-resolution-type finite difference (or finite volume) method of the family introduced by Wang, Shu, Yee and Sjögreen (2012) and extended to multispecies, multireaction reactive Euler flows by Wang, Shu, Yee, Kotov and Sjögreen (2015), for advection–reaction equations $u_t + f(u)_x = s(u)$ with stiff reaction terms — and for reactive Euler systems with the multi-step, multi-species kinetics classes treated in those SR papers — that, on meshes which grossly under-resolve the reaction zone (one or a few grid points inside it), simultaneously achieves: (1) correct-speed discontinuity capturing, i.e., under-resolved reaction fronts appear at the correct locations with the correct jump sizes (correct propagation speeds), the capability the 2012 and 2015 subcell-resolution methods already provide and must retain; (2) high order in time without operator splitting, i.e., at least third-order temporal accuracy, removing the second-order cap that Strang splitting places on the existing SR method, with the troubled-point subcell correction embedded consistently in every Runge–Kutta stage or multistep level of the non-splitting integrator; and (3) implicit compatibility, i.e., the stiff reaction step can be integrated by implicit or IMEX schemes, so the troubled-point correction must no longer preclude implicit treatment of the stiff dynamics — in the existing SR method, implicit methods cannot be used in the reaction step because the troubled values need to be modified explicitly. Capabilities (2) and (3) correspond to limitations (2) and (4) acknowledged by the authors in their condensed CTR 2011 text. The construction defines the algorithm (troubled-cell detection, subcell reconstruction, time discretization), and its accuracy on coarse-mesh stiff advection–reaction and reactive-Euler detonation problems is established against fine-mesh resolved reference solutions. The multi-species, multi-reaction extension of subcell resolution is not itself novel — it was resolved by the same authors in 2015, albeit still within a fractional-step framework with an explicitly modified troubled-value reaction step — but capabilities (2) and (3) must also be demonstrated for those multispecies, multireaction kinetics. Resolution of the internal structure of the reaction zone is likewise outside the coarse-mesh objective: as the authors remark, resolving the sharp reaction zone requires sufficiently many grid points in that zone.

References

1. Wang, Shu, Yee, Sjögreen, 'High order finite difference methods with subcell resolution for advection equations with stiff source terms', Journal of Computational Physics 231(1):190–214, 2012. DOI 10.1016/j.jcp.2011.08.031, https://doi.org/10.1016/j.jcp.2011.08.031
2. Wang, Shu, Yee, Kotov, Sjögreen, 'High Order Finite Difference Methods with Subcell Resolution for Stiff Multispecies Discontinuity Capturing', Communications in Computational Physics 17(2):317–336, 2015. DOI 10.4208/cicp.250214.130814a, https://doi.org/10.4208/cicp.250214.130814a
3. Wang, Shu, Yee, Sjögreen, 'High order finite difference methods with subcell resolution for 2D detonation waves', Annual Research Briefs 2011, Center for Turbulence Research, Stanford University. NTRS accession 20120011106, https://ntrs.nasa.gov/citations/20120011106
4. Du, Wang, Qian, Yang, 'High-Order Bound-Preserving Discontinuous Galerkin Methods for Stiff Multispecies Detonation', SIAM Journal on Scientific Computing 41(2):B250–B273, 2019. DOI 10.1137/18m122265x, https://doi.org/10.1137/18m122265x
5. Du, Yang, 'Third-order conservative sign-preserving and steady-state-preserving time integrations and applications in stiff multispecies and multireaction detonations', Journal of Computational Physics 395:489–510, 2019. DOI 10.1016/j.jcp.2019.06.040, https://doi.org/10.1016/j.jcp.2019.06.040
6. Du, Yang, 'High-Order Bound-Preserving Finite Difference Methods for Multispecies and Multireaction Detonations', Communications on Applied Mathematics and Computation 5:31–63, 2021. DOI 10.1007/s42967-020-00117-y, https://doi.org/10.1007/s42967-020-00117-y
7. Huang, Shu, 'Bound-preserving modified exponential Runge–Kutta discontinuous Galerkin methods for scalar hyperbolic equations with stiff source terms', Journal of Computational Physics 361:111–135, 2018. DOI 10.1016/j.jcp.2018.01.051, https://doi.org/10.1016/j.jcp.2018.01.051
8. Chertock, Chu, Kurganov, 'Accurate deterministic projection methods for stiff detonation waves', Communications in Mathematical Sciences 22(4):871–910, 2024. DOI 10.4310/cms.2024.v22.n4.a1, https://doi.org/10.4310/cms.2024.v22.n4.a1
9. Wang, Pan, Hu, Adams, 'A split random time-stepping method for stiff and nonstiff detonation capturing', Combustion and Flame 204:397–413, 2019 (preprint arXiv:1802.04116, titled 'A Split Random Time Stepping Method for Stiff and Non-stiff Chemically Reacting Flows'). DOI 10.1016/j.combustflame.2019.03.034, https://doi.org/10.1016/j.combustflame.2019.03.034
10. Deng, Xie, Loubère, Shimizu, Xiao, 'Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts', Computers and Fluids 171:1–14, 2018. DOI 10.1016/j.compfluid.2018.05.015, https://doi.org/10.1016/j.compfluid.2018.05.015
11. Deng, Xie, Xiao, Teng, 'New Accurate and Efficient Method for Stiff Detonation Capturing', AIAA Journal 56(10):4024–4038, 2018. DOI 10.2514/1.J056632, https://doi.org/10.2514/1.J056632
12. Yee, Kotov, Wang, Shu, 'Spurious behavior of shock-capturing methods by the fractional step approach: Problems containing stiff source terms and discontinuities', Journal of Computational Physics 241:266–291, 2013. DOI 10.1016/j.jcp.2013.01.028, https://doi.org/10.1016/j.jcp.2013.01.028
13. Zhu, Huang, Yang, 'Bound-preserving discontinuous Galerkin methods with modified Patankar time integrations for chemical reacting flows', preprint, 2022. arXiv:2211.05600, https://arxiv.org/abs/2211.05600
14. Huo, 'PDE Realization and Structure-Aware Solvers for a Compact Two-Stage Fourth-Order IMEX Method', preprint, 2026 (sole author Zhixin Huo). arXiv:2608.13875, https://arxiv.org/abs/2608.13875