On Arbitrary-Lagrangian-Eulerian One-Step WENO Schemes for Stiff Hyperbolic Balance Laws

In this article we present a new family of high order accurate Arbitrary Lagrangian-Eulerian one-step WENO finite volume schemes for the solution of stiff hyperbolic balance laws. High order accuracy in space is obtained with a standard WENO reconstruction algorithm and high order in time is obtained using the local space-time discontinuous Galerkin method recently proposed in [20]. In the Lagrangian framework considered here, the local space-time DG predictor is based on a weak formulation of the governing PDE on a moving space-time element. For the space-time basis and test functions we use Lagrange interpolation polynomials defined by tensor-product Gauss-Legendre quadrature points. The moving space-time elements are mapped to a reference element using an isoparametric approach, i.e. the space-time mapping is defined by the same basis functions as the weak solution of the PDE. We show some computational examples in one space-dimension for non-stiff and for stiff balance laws, in particular for the Euler equations of compressible gas dynamics, for the resistive relativistic MHD equations, and for the relativistic radiation hydrodynamics equations. Numerical convergence results are presented for the stiff case up to sixth order of accuracy in space and time and for the non-stiff case up to eighth order of accuracy in space and time.     

Arbitrary-Lagrangian-Eulerian One-Step WENO Finite Volume Schemes on Unstructured Triangular Meshes

In this article we present a new class of high order accurate ArbitraryEulerian-Lagrangian (ALE) one-step WENO finite volume schemes for solving nonlinear hyperbolic systems of conservation laws on moving two dimensional unstructured triangular meshes. A WENO reconstruction algorithm is used to achieve high order accuracy in space and a high order one-step time discretization is achieved by using the local space-time Galerkin predictor proposed in [25]. For that purpose, a new element-local weak formulation of the governing PDE is adopted on moving space-time elements. The space-time basis and test functions are obtained considering Lagrange interpolation polynomials passing through a predefined set of nodes. Moreover, a polynomial mapping defined by the same local space-time basis functions as the weak solution of the PDE is used to map the moving physical space-time element onto a space-time reference element. To maintain algorithmic simplicity, the final ALE one-step finite volume scheme uses moving triangular meshes with straight edges. This is possible in the ALE framework, which allows a local mesh velocity that is different from the local fluid velocity. We present numerical convergence rates for the schemes presented in this paper up to sixth order of accuracy in space and time and show some classical numerical test problems for the two-dimensional Euler equations of compressible gas dynamics.     

Numerical Methods for Balance Laws with Space Dependent Flux: Application to Radiotherapy Dose Calculation

The present work is concerned with the derivation of numerical methods to approximate the radiation dose in external beam radiotherapy. To address this issue, we consider a moment approximation of radiative transfer, closed by an entropy minimization principle. The model under consideration is governed by a system of hyperbolic equations in conservation form supplemented by source terms. The main difficulty coming from the numerical approximation of this system is an explicit space dependence in the flux function. Indeed, this dependence will be seen to be stiff and specific numerical strategies must be derived in order to obtain the needed accuracy. A first approach is developed considering the 1D case, where a judicious change of variables allows to eliminate the space dependence in the flux function. This is not possible in multi-D. We therefore reinterpret the 1D scheme as a scheme on two meshes, and generalize this to 2D by alternating transformations between separate meshes. We call this procedure projection method. Several numerical experiments, coming from medical physics, illustrate the potential applicability of the developed method.     

A Compact High Order Space-Time Method for Conservation Laws

This paper presents a novel high-order space-time method for hyperbolic conservation laws. Two important concepts, the staggered space-time mesh of the space-time conservation element/solution element (CE/SE) method and the local discontinuous basis functions of the space-time discontinuous Galerkin (DG) finite element method, are the two key ingredients of the new scheme. The staggered space-time mesh is constructed using the cell-vertex structure of the underlying spatial mesh. The universal definitions of CEs and SEs are independent of the underlying spatial mesh and thus suitable for arbitrarily unstructured meshes. The solution within each physical time step is updated alternately at the cell level and the vertex level. For this solution updating strategy and the DG ingredient, the new scheme here is termed as the discontinuous Galerkin cell-vertex scheme (DG-CVS). The high order of accuracy is achieved by employing high-order Taylor polynomials as the basis functions inside each SE. The present DG-CVS exhibits many advantageous features such as Riemann-solver-free, high-order accuracy, point-implicitness, compactness, and ease of handling boundary conditions. Several numerical tests including the scalar advection equations and compressible Euler equations will demonstrate the performance of the new method.     

An All-Speed Asymptotic-Preserving Method for the Isentropic Euler and Navier-Stokes Equations

The computation of compressible flows becomes more challenging when the Mach number has different orders of magnitude. When the Mach number is of order one, modern shock capturing methods are able to capture shocks and other complex structures with high numerical resolutions. However, if the Mach number is small, the acoustic waves lead to stiffness in time and excessively large numerical viscosity, thus demanding much smaller time step and mesh size than normally needed for incompressible flow simulation. In this paper, we develop an all-speed asymptotic preserving (AP) numerical scheme for the compressible isentropic Euler and Navier-Stokes equations that is uniformly stable and accurate for all Mach numbers. Our idea is to split the system into two parts: one involves a slow, nonlinear and conservative hyperbolic system adequate for the use of modern shock capturing methods and the other a linear hyperbolic system which contains the stiff acoustic dynamics, to be solved implicitly. This implicit part is reformulated into a standard pressure Poisson projection system and thus possesses sufficient structure for efficient fast Fourier transform solution techniques. In the zero Mach number limit, the scheme automatically becomes a projection method-like incompressible solver. We present numerical results in one and two dimensions in both compressible and incompressible regimes.     

Continuous Finite Element Subgrid Basis Functions for Discontinuous Galerkin Schemes on Unstructured Polygonal Voronoi Meshes

We propose a new high order accurate nodal discontinuous Galerkin (DG) method for the solution of nonlinear hyperbolic systems of partial differential equations (PDE) on unstructured polygonal Voronoi meshes. Rather than using classical polynomials of degree $N$ inside each element, in our new approach the discrete solution is represented by piecewise continuous polynomials of degree $N$ within each Voronoi element, using a continuous finite element basis defined on a subgrid inside each polygon. We call the resulting subgrid basis an agglomerated finite element (AFE) basis for the DG method on general polygons, since it is obtained by the agglomeration of the finite element basis functions associated with the subgrid triangles. The basis functions on each sub-triangle are defined, as usual, on a universal reference element, hence allowing to compute universal mass, flux and stiffness matrices for the subgrid triangles once and for all in a pre-processing stage for the reference element only. Consequently, the construction of an efficient quadrature-free algorithm is possible, despite the unstructured nature of the computational grid. High order of accuracy in time is achieved thanks to the ADER approach, making use of an element-local space-time Galerkin finite element predictor.The novel schemes are carefully validated against a set of typical benchmark problems for the compressible Euler and Navier-Stokes equations. The numerical results have been checked with reference solutions available in literature and also systematically compared, in terms of computational efficiency and accuracy, with those obtained by the corresponding modal DG version of the scheme.     

A New Approach of High Order Well-Balanced Finite Volume WENO Schemes and Discontinuous Galerkin Methods for a Class of Hyperbolic Systems with Source Terms

Hyperbolic balance laws have steady state solutions in which the flux gradients are nonzero but are exactly balanced by the source terms. In our earlier work [31–33], we designed high order well-balanced schemes to a class of hyperbolic systems with separable source terms. In this paper, we present a different approach to the same purpose: designing high order well-balanced finite volume weighted essentially non-oscillatory (WENO) schemes and RungeKutta discontinuous Galerkin (RKDG) finite element methods. We make the observation that the traditional RKDG methods are capable of maintaining certain steady states exactly, if a small modification on either the initial condition or the flux is provided. The computational cost to obtain such a well balanced RKDG method is basically the same as the traditional RKDG method. The same idea can be applied to the finite volume WENO schemes. We will first describe the algorithms and prove the well balanced property for the shallow water equations, and then show that the result can be generalized to a class of other balance laws. We perform extensive one and two dimensional simulations to verify the properties of these schemes such as the exact preservation of the balance laws for certain steady state solutions, the non-oscillatory property for general solutions with discontinuities, and the genuine high order accuracy in smooth regions.     

On Universal Osher-Type Schemes for General Nonlinear Hyperbolic Conservation Laws

This paper is concerned with a new version of the Osher-Solomon Riemann solver and is based on a numerical integration of the path-dependent dissipation matrix. The resulting scheme is much simpler than the original one and is applicable to general hyperbolic conservation laws, while retaining the attractive features of the original solver: the method is entropy-satisfying, differentiable and complete in the sense that it attributes a different numerical viscosity to each characteristic field, in particular to the intermediate ones, since the full eigenstructure of the underlying hyperbolic system is used. To illustrate the potential of the proposed scheme we show applications to the following hyperbolic conservation laws: Euler equations of compressible gasdynamics with ideal gas and real gas equation of state, classical and relativistic MHD equations as well as the equations of nonlinear elasticity. To the knowledge of the authors, apart from the Euler equations with ideal gas, an Osher-type scheme has never been devised before for any of these complicated PDE systems. Since our new general Riemann solver can be directly used as a building block of high order finite volume and discontinuous Galerkin schemes we also show the extension to higher order of accuracy and multiple space dimensions in the new framework of P_NP_M schemes on unstructured meshes recently proposed in [9].     

Relaxation methods for hyperbolic PDE mixed‐integer optimal control problems

The convergence analysis for methods solving partial differential equations constrained optimal control problems containing both discrete and continuous control decisions based on relaxation and rounding strategies is extended to the class of first order semilinear hyperbolic systems in one space dimension. The results are obtained by novel a priori estimates for the size of the relaxation gap based on the characteristic flow, fixed‐point arguments, and particular regularity theory for such mixed‐integer control problems. Motivated by traffic flow problems, a relaxation model for optimal flux switching control in conservation laws is considered as an application.     

Extension of the High-Order Space-Time Discontinuous Galerkin Cell Vertex Scheme to Solve Time Dependent Diffusion Equations

In this paper, the high-order space-time discontinuous Galerkin cell vertex scheme (DG-CVS) developed by the authors for hyperbolic conservation laws is extended for time dependent diffusion equations. In the extension, the treatment of the diffusive flux is exactly the same as that for the advective flux. Thanks to the Riemann-solver-free and reconstruction-free features of DG-CVS, both the advective flux and the diffusive flux are evaluated using continuous information across the cell interface. As a result, the resulting formulation with diffusive fluxes present is still consistent and does not need any extra ad hoc techniques to cure the common "variational crime" problem when traditional DG methods are applied to diffusion problems. For this reason, DG-CVS is conceptually simpler than other existing DG-typed methods. The numerical tests demonstrate that the convergence order based on the L²-norm is optimal, i.e. O(h^p+1) for the solution and O(h^p) for the solution gradients, when the basis polynomials are of odd degrees. For even-degree polynomials, the convergence order is sub-optimal for the solution and optimal for the solution gradients. The same odd-even behaviour can also be seen in some other DG-typed methods.     

Relaxation Schemes for the $M_1$ Model with Space-Dependent Flux: Application to Radiotherapy Dose Calculation

Because of stability constraints, most numerical schemes applied to hyperbolic systems of equations turn out to be costly when the flux term is multiplied by some very large scalar. This problem emerges with the $M_1$ system of equations in the field of radiotherapy when considering heterogeneous media with very disparate densities. Additionally, the flux term of the $M_1$ system is non-linear, and in order for the model to be well-posed the numerical solution needs to fulfill conditions called realizability. In this paper, we propose a numerical method that overcomes the stability constraint and preserves the realizability property. For this purpose, we relax the $M_1$ system to obtain a linear flux term. Then we extend the stencil of the difference quotient to obtain stability. The scheme is applied to a radiotherapy dose calculation example.     

A Well-Balanced Gas Kinetic Scheme for Navier-Stokes Equations with Gravitational Potential

The hydrostatic equilibrium state is the consequence of the exact balance between hydrostatic pressure and external force. Standard finite volume cannot keep this balance exactly due to their unbalanced truncation errors. In this study, we introduce an auxiliary variable which becomes constant at isothermal hydrostatic equilibria and propose a well-balanced gas kinetic scheme for the Navier-Stokes equations. Through reformulating the convection term and the force term via the auxiliary variable, zero numerical flux and zero numerical source term are enforced at the hydrostatic equilibrium state instead of the balance between hydrostatic pressure and external force. Several problems are tested to demonstrate the accuracy and the stability of the new scheme. The results confirm that, the new scheme can preserve the exact hydrostatic solution. The small perturbation riding on hydrostatic equilibria can be calculated accurately. More importantly, the new scheme is capable of simulating the process of converging towards hydrostatic equilibria from a highly unbalanced initial condition. The ultimate state of zero velocity and constant temperature is achieved up to machine accuracy. As demonstrated by the numerical experiments, the current scheme is very suitable for small amplitude perturbation and long time running under gravitational potential.     

Numerical Entropy and Adaptivity for Finite Volume Schemes

We propose an a-posteriori error/smoothness indicator for standard semi-discrete finite volume schemes for systems of conservation laws, based on the numerical production of entropy. This idea extends previous work by the first author limited to central finite volume schemes on staggered grids. We prove that the indicator converges to zero with the same rate of the error of the underlying numerical scheme on smooth flows under grid refinement. We construct and test an adaptive scheme for systems of equations in which the mesh is driven by the entropy indicator. The adaptive scheme uses a single nonuniform grid with a variable timestep. We show how to implement a second order scheme on such a space-time non uniform grid, preserving accuracy and conservation properties. We also give an example of a p-adaptive strategy.     

A High Order Sharp-Interface Method with Local Time Stepping for Compressible Multiphase Flows

In this paper, a new sharp-interface approach to simulate compressible multiphase flows is proposed. The new scheme consists of a high order WENO finite volume scheme for solving the Euler equations coupled with a high order path-conservative discontinuous Galerkin finite element scheme to evolve an indicator function that tracks the material interface. At the interface our method applies ghost cells to compute the numerical flux, as the ghost fluid method. However, unlike the original ghost fluid scheme of Fedkiw et al. [15], the state of the ghost fluid is derived from an approximate-state Riemann solver, similar to the approach proposed in [25], but based on a much simpler formulation. Our formulation leads only to one single scalar nonlinear algebraic equation that has to be solved at the interface, instead of the system used in [25]. Away from the interface, we use the new general Osher-type flux recently proposed by Dumbser and Toro [13], which is a simple but complete Riemann solver, applicable to general hyperbolic conservation laws. The time integration is performed using a fully-discrete one-step scheme, based on the approaches recently proposed in [5, 7]. This allows us to evolve the system also with time-accurate local time stepping. Due to the sub-cell resolution and the subsequent more restrictive time-step constraint of the DG scheme, a local evolution for the indicator function is applied, which is matched with the finite volume scheme for the solution of the Euler equations that runs with a larger time step. The use of a locally optimal time step avoids the introduction of excessive numerical diffusion in the finite volume scheme. Two different fluids have been used, namely an ideal gas and a weakly compressible fluid modeled by the Tait equation. Several tests have been computed to assess the accuracy and the performance of the new high order scheme. A verification of our algorithm has been carefully carried out using exact solutions as well as a comparison with other numerical reference solutions. The material interface is resolved sharply and accurately without spurious oscillations in the pressure field.     

A Lattice Boltzmann and Immersed Boundary Scheme for Model Blood Flow in Constricted Pipes: Part 1 – Steady Flow

Hemodynamics is a complex problem with several distinct characteristics; fluid is non-Newtonian, flow is pulsatile in nature, flow is three-dimensional due to cholesterol/plague built up, and blood vessel wall is elastic. In order to simulate this type of flows accurately, any proposed numerical scheme has to be able to replicate these characteristics correctly, efficiently, as well as individually and collectively. Since the equations of the finite difference lattice Boltzmann method (FDLBM) are hyperbolic, and can be solved using Cartesian grids locally, explicitly and efficiently on parallel computers, a program of study to develop a viable FDLBM numerical scheme that can mimic these characteristics individually in any model blood flow problem was initiated. The present objective is to first develop a steady FDLBM with an immersed boundary (IB) method to model blood flow in stenoic artery over a range of Reynolds numbers. The resulting equations in the FDLBM/IB numerical scheme can still be solved using Cartesian grids; thus, changing complex artery geometry can be treated without resorting to grid generation. The FDLBM/IB numerical scheme is validated against known data and is then used to study Newtonian and non-Newtonian fluid flow through constricted tubes. The investigation aims to gain insight into the constricted flow behavior and the non-Newtonian fluid effect on this behavior.     

A Compact Third-Order Gas-Kinetic Scheme for Compressible Euler and Navier-Stokes Equations

In this paper, a compact third-order gas-kinetic scheme is proposed for the compressible Euler and Navier-Stokes equations. The main reason for the feasibility to develop such a high-order scheme with compact stencil, which involves only neighboring cells, is due to the use of a high-order gas evolution model. Besides the evaluation of the time-dependent flux function across a cell interface, the high-order gas evolution model also provides an accurate time-dependent solution of the flow variables at a cell interface. Therefore, the current scheme not only updates the cell averaged conservative flow variables inside each control volume, but also tracks the flow variables at the cell interface at the next time level. As a result, with both cell averaged and cell interface values, the high-order reconstruction in the current scheme can be done compactly. Different from using a weak formulation for high-order accuracy in the Discontinuous Galerkin method, the current scheme is based on the strong solution, where the flow evolution starting from a piecewise discontinuous high-order initial data is precisely followed. The cell interface time-dependent flow variables can be used for the initial data reconstruction at the beginning of next time step. Even with compact stencil, the current scheme has third-order accuracy in the smooth flow regions, and has favorable shock capturing property in the discontinuous regions. It can be faithfully used from the incompressible limit to the hypersonic flow computations, and many test cases are used to validate the current scheme. In comparison with many other high-order schemes, the current method avoids the use of Gaussian points for the flux evaluation along the cell interface and the multi-stage Runge-Kutta time stepping technique. Due to its multidimensional property of including both derivatives of flow variables in the normal and tangential directions of a cell interface, the viscous flow solution, especially those with vortex structure, can be accurately captured. With the same stencil of a second order scheme, numerical tests demonstrate that the current scheme is as robust as well-developed second-order shock capturing schemes, but provides more accurate numerical solutions than the second order counterparts.     

A One-Dimensional Second-Order Cell-Centered Lagrangian Scheme Satisfying the Entropy Condition

The numerical solutions of gas dynamics equations have to be consistent with the second law of thermodynamics, which is termed entropy condition. However, most cell-centered Lagrangian (CL) schemes do not satisfy the entropy condition. Until 2020, for one-dimensional gas dynamics equations, the first-order CL scheme with the hybridized flux developed by combining the acoustic approximate (AA) flux and the entropy conservative (EC) flux developed by Maire et al. was used. This hybridized CL scheme satisfies the entropy condition; however, it is under-entropic in the part zones of rarefaction waves. Moreover, the EC flux may result in nonphysical numerical oscillations in simulating strong rarefaction waves. Another disadvantage of this scheme is that it is of only first-order accuracy. In this paper, we firstly construct a modified entropy conservative (MEC) flux which can damp effectively numerical oscillations in simulating strong rarefaction waves. Then we design a new hybridized CL scheme satisfying the entropy condition for one-dimensional complex flows. This new hybridized CL scheme is a combination of the AA flux and the MEC flux.In order to prevent the specific entropy of the hybridized CL scheme from being under-entropic, we propose using the third-order TVD-type Runge-Kutta time discretization method. Based on the new hybridized flux, we develop the second-order CL scheme that satisfies the entropy condition.Finally, the characteristics of our new CL scheme using the improved hybridized flux are demonstrated through several numerical examples.     

A New Family of High Order Unstructured MOOD and ADER Finite Volume Schemes for Multidimensional Systems of Hyperbolic Conservation Laws

In this paper, we investigate the coupling of the Multi-dimensional Optimal Order Detection (MOOD) method and the Arbitrary high order DERivatives (ADER) approach in order to design a new high order accurate, robust and computationally efficient Finite Volume (FV) scheme dedicated to solving nonlinear systems of hyperbolic conservation laws on unstructured triangular and tetrahedral meshes in two and three space dimensions, respectively. The Multi-dimensional Optimal Order Detection (MOOD) method for 2D and 3D geometries has been introduced in a recent series of papers for mixed unstructured meshes. It is an arbitrary high-order accurate Finite Volume scheme in space, using polynomial reconstructions with a posteriori detection and polynomial degree decrementing processes to deal with shock waves and other discontinuities. In the following work, the time discretization is performed with an elegant and efficient one-step ADER procedure. Doing so, we retain the good properties of the MOOD scheme, that is to say, the optimal high-order of accuracy is reached on smooth solutions, while spurious oscillations near singularities are prevented. The ADER technique not only reduces the cost of the overall scheme as shown on a set of numerical tests in 2D and 3D, but also increases the stability of the overall scheme. A systematic comparison between classical unstructured ADER-WENO schemes and the new ADER-MOOD approach has been carried out for high-order schemes in space and time in terms of cost, robustness, accuracy and efficiency. The main finding of this paper is that the combination of ADER with MOOD generally outperforms the one of ADER and WENO either because at given accuracy MOOD isless expensive (memory and/or CPU time), or because it is more accurate for a given grid resolution. A large suite of classical numerical test problems has been solved on unstructured meshes for three challenging multi-dimensional systems of conservation laws: the Euler equations of compressible gas dynamics, the classical equations of ideal magneto-Hydrodynamics (MHD) and finally the relativistic MHD equations (RMHD), which constitutes a particularly challenging nonlinear system of hyperbolic partial differential equation. All tests are run on genuinely unstructured grids composed of simplex elements.     

Immersed Interface CIP for One Dimensional Hyperbolic Equations

The immersed interface technique is incorporated into CIP method to solve one-dimensional hyperbolic equations with piecewise constant coefficients. The proposed method achieves the third order of accuracy in time and space in the vicinity of the interface where the coefficients have jump discontinuities, which is the same order of accuracy of the standard CIP scheme. Some numerical tests are given to verify the accuracy of the proposed method.     

All Speed Scheme for the Low Mach Number Limit of the Isentropic Euler Equations

An all speed scheme for the Isentropic Euler equations is presented in this paper. When the Mach number tends to zero, the compressible Euler equations converge to their incompressible counterpart, in which the density becomes a constant. Increasing approximation errors and severe stability constraints are the main difficulty in the low Mach regime. The key idea of our all speed scheme is the special semi-implicit time discretization, in which the low Mach number stiff term is divided into two parts, one being treated explicitly and the other one implicitly. Moreover, the flux of the density equation is also treated implicitly and an elliptic type equation is derived to obtain the density. In this way, the correct limit can be captured without requesting the mesh size and time step to be smaller than the Mach number. Compared with previous semi-implicit methods [11,13,29], firstly, nonphysical oscillations can be suppressed by choosing proper parameter, besides, only a linear elliptic equation needs to be solved implicitly which reduces much computational cost. We develop this semi-implicit time discretization in the framework of a first order Local Lax-Friedrichs (or Rusanov) scheme and numerical tests are displayed to demonstrate its performances.