In this paper, the author recounts his forty-year plus struggle to find a sound basis for understanding the computational fluid dynamics of compressible flow.

In this paper, we discuss the notion of discrete conservation for hyperbolic conservation laws. We introduce what we call fluctuation splitting schemes (or residual distribution, also RDS) and show through several examples how these schemes lead to new developments. In particular, we show that most, if not all, known schemes can be rephrased in flux form and also how to satisfy additional conservation laws. This review paper is built on Abgrall et al. (Computers and Fluids 169:10–22, 2018), Abgrall and Tokareva (SIAM SISC 39(5):A2345–A2364, 2017), Abgrall (J Sci Comput 73:461–494, 2017), Abgrall (Methods Appl Math 18(3):327–351, 2018a) and Abgrall (J Comput Phys 372, 640–666, 2018b). This paper is also a direct consequence of the work of Roe, in particular Deconinck et al. (Comput Fluids 22(2/3):215–222, 1993) and Roe (J Comput Phys 43:357–372, 1981) where the notion of conservation was first introduced. In [26], Roe mentioned the Hermes project and the role of Dassault Aviation. Bruno Stoufflet, Vice President R&D and advanced business of this company, proposed me to have a detailed look at Deconinck et al. (Comput Fluids 22(2/3):215–222, 1993). To be honest, at the time, I did not understand anything, and this was the case for several years. I was lucky to work with Katherine Mer, who at the time was a postdoc, and is now research engineer at CEA. She helped me a lot in understanding the notion of conservation. The present contribution can be seen as the result of my understanding after many years of playing around with the notion of residual distribution schemes (or fluctuation-splitting schemes) introduced by Roe.

The ADER approach to solve hyperbolic equations to very high order of accuracy has seen explosive developments in the last few years, including both methodological aspects as well as very ambitious applications. In spite of methodological progress, the issues of efficiency and ease of implementation of the solution of the associated generalized Riemann problem (GRP) remain the centre of attention in the ADER approach. In the original formulation of ADER schemes, the proposed solution procedure for the GRP was based on (i) Taylor series expansion of the solution in time right at the element interface, (ii) subsequent application of the Cauchy–Kowalewskaya procedure to convert time derivatives to functionals of space derivatives, and (iii) solution of classical Riemann problems for high-order spatial derivatives to complete the Taylor series expansion. For realistic problems the Cauchy–Kowalewskaya procedure requires the use of symbolic manipulators and being rather cumbersome its replacement or simplification is highly desirable. In this paper we propose a new class of solvers for the GRP that avoid the Cauchy–Kowalewskaya procedure and result in simpler ADER schemes. This is achieved by exploiting the history of the numerical solution that makes it possible to devise a time-reconstruction procedure at the element interface. Still relying on a time Taylor series expansion of the solution at the interface, the time derivatives are then easily calculated from the time-reconstruction polynomial. The resulting schemes are called ADER-TR. A thorough study of the linear stability properties of the linear version of the schemes is carried out using the von Neumann method, thus deducing linear stability regions. Also, via careful numerical experiments, we deduce stability regions for the corresponding non-linear schemes. Numerical examples using the present simplified schemes of fifth and seventh order of accuracy in space and time show that these compare favourably with conventional ADER methods. This paper is restricted to the one-dimensional scalar case with source term, but preliminary results for the one-dimensional Euler equations indicate that the time-reconstruction approach offers significant advantages not only in terms of ease of implementation but also in terms of efficiency for the high-order range schemes.

A new type of high-order multi-resolution weighted essentially non-oscillatory (WENO) schemes (Zhu and Shu in J Comput Phys, 375: 659–683, 2018) is applied to solve for steady-state problems on structured meshes. Since the classical WENO schemes (Jiang and Shu in J Comput Phys, 126: 202–228, 1996) might suffer from slight post-shock oscillations (which are responsible for the residue to hang at a truncation error level), this new type of high-order finite-difference and finite-volume multi-resolution WENO schemes is applied to control the slight post-shock oscillations and push the residue to settle down to machine zero in steady-state simulations. This new type of multi-resolution WENO schemes uses the same large stencils as that of the same order classical WENO schemes, could obtain fifth-order, seventh-order, and ninth-order in smooth regions, and could gradually degrade to first-order so as to suppress spurious oscillations near strong discontinuities. The linear weights of such new multi-resolution WENO schemes can be any positive numbers on the condition that their sum is one. This is the first time that a series of unequal-sized hierarchical central spatial stencils are used in designing high-order finite-difference and finite-volume WENO schemes for solving steady-state problems. In comparison with the classical fifth-order finite-difference and finite-volume WENO schemes, the residue of these new high-order multi-resolution WENO schemes can converge to a tiny number close to machine zero for some benchmark steady-state problems.

One of the beneficial properties of the discontinuous Galerkin method is the accurate wave propagation properties. That is, the semi-discrete error has dissipation errors of order $2k+1$ ($\le Ch^{2k+1}$) and order $2k+2$ for dispersion ($\le Ch^{2k+2}$). Previous studies have concentrated on the order of accuracy, and neglected the important role that the error constant, C,  plays in these estimates. In this article, we show the important role of the error constant in the dispersion and dissipation error for discontinuous Galerkin approximation of polynomial degree k,  where $k=0,1,2,3.$ This gives insight into why one may want a more centred flux for a piecewise constant or quadratic approximation than for a piecewise linear or cubic approximation. We provide an explicit formula for these error constants. This is illustrated through one particular flux, the upwind-biased flux introduced by Meng et al., as it is a convex combination of the upwind and downwind fluxes. The studies of wave propagation are typically done through a Fourier ansatz. This higher order Fourier information can be extracted using the smoothness-increasing accuracy-conserving (SIAC) filter. The SIAC filter ties the higher order Fourier information to the negative-order norm in physical space. We show that both the proofs of the ability of the SIAC filter to extract extra accuracy and numerical results are unaffected by the choice of flux.

We propose a numerical methodology for the simultaneous numerical simulation of four states of matter: gas, liquid, elastoplastic solids, and plasma. The distinct, interacting physical processes are described by a combination of compressible, inert, and reactive forms of the Euler equations, multi-phase equations, elastoplastic equations, and resistive MHD equations. Combinations of systems of equations are usually solved by coupling finite element for solid modelling and CFD models for fluid modelling or including material effects through boundary conditions rather than full material discretisation. Our simultaneous solution methodology lies on the recasting of all the equations in the same, hyperbolic form allowing their solution on the same grid with the same finite volume numerical schemes. We use a combination of sharp- and diffuse-interface methods to track or capture material interfaces, depending on the application. The communication between the distinct systems of equations (i.e., materials separated by sharp interfaces) is facilitated by means of mixed-material Riemann solvers at the boundaries of the systems, which represent physical material boundaries. To this end, we derive approximate mixed-material Riemann solvers for each pair of the above models based on characteristic equations. To demonstrate the applicability of the new methodology, we consider a case study, where we investigate the possibility of ignition of a combustible gas that lies over a liquid in a metal container that is struck by a plasma arc akin to a lightning strike. We study the effect of the metal container material and its conductivity on the ignition of the combustible gas, as well as the effects of an additional dielectric coating, the sensitivity of the gas, and differences between scenarios with sealed and pre-damaged metal surfaces.

We construct an approximate Riemann solver for scalar advection–diffusion equations with piecewise polynomial initial data. The objective is to handle advection and diffusion simultaneously to reduce the inherent numerical diffusion produced by the usual advection flux calculations. The approximate solution is based on the weak formulation of the Riemann problem and is solved within a space–time discontinuous Galerkin approach with two subregions. The novel generalized Riemann solver produces piecewise polynomial solutions of the Riemann problem. In conjunction with a recovery polynomial, the Riemann solver is then applied to define the numerical flux within a finite volume method. Numerical results for a piecewise linear and a piecewise parabolic approximation are shown. These results indicate a reduction in numerical dissipation compared with the conventional separated flux calculation of advection and diffusion. Also, it is shown that using the proposed solver only in the vicinity of discontinuities gives way to an accurate and efficient finite volume scheme.