James A. Rossmanith

Camille Felton, Mariana Harris, Caleb Logemann et al.

Nonlinear hyperbolic conservation laws admit singular solutions such as shockwaves (discontinuities in conserved variables), rarefaction waves (discontinuities in derivatives), and vacuum states (loss of strong hyperbolicity). When ostensibly high-order numerical methods are applied in such solution regimes, unphysical oscillations present themselves that can lead to large errors and a breakdown of the numerical simulation. In this work we develop a new Lax-Wendroff discontinuous Galerkin (LxW-DG) method with a limiting strategy that keeps the solution non-oscillatory and positivity-preserving for relevant variables, such as height in the shallow water equations and density and pressure in the compressible Euler equations. The proposed LxW-DG scheme updates the solution over each time-step with a locally-implicit predictor followed by an explicit corrector. The locally-implicit prediction phase is formulated in terms of primitive variables, which greatly simplifies the solver. The resulting system of nonlinear algebraic equations are approximately solved via a Picard iteration, where the number of iterations is equal to the order of accuracy of the method. The correction phase is an explicit evaluation formulated in terms of conservative variables in order to guarantee numerical conservation. In order to achieve full positivity-preservation, limiting is required in both the prediction and correction steps. The resulting scheme is applied to several standard test cases for the shallow water and compressible Euler equations. All of the presented examples are written in a freely available open-source Python code.

Christiane Helzel, James A. Rossmanith, Bertram Taetz

Numerical methods for solving the ideal magnetohydrodynamic (MHD) equations in more than one space dimension must either confront the challenge of controlling errors in the discrete divergence of the magnetic field, or else be faced with nonlinear numerical instabilities. One approach for controlling the discrete divergence is through a so-called constrained transport method, which is based on first predicting a magnetic field through a standard finite volume solver, and then correcting this field through the appropriate use of a magnetic vector potential. In this work we develop a constrained transport method for the 3D ideal MHD equations that is based on a high-resolution wave propagation scheme. Our proposed scheme is the 3D extension of the 2D scheme developed by Rossmanith [SIAM J. Sci. Comp. 28, 1766 (2006)], and is based on the high-resolution wave propagation method of Langseth and LeVeque [J. Comp. Phys. 165, 126 (2000)]. In particular, in our extension we take great care to maintain the three most important properties of the 2D scheme: (1) all quantities, including all components of the magnetic field and magnetic potential, are treated as cell-centered; (2) we develop a high-resolution wave propagation scheme for evolving the magnetic potential; and (3) we develop a wave limiting approach that is applied during the vector potential evolution, which controls unphysical oscillations in the magnetic field. One of the key numerical difficulties that is novel to 3D is that the transport equation that must be solved for the magnetic vector potential is only weakly hyperbolic. In presenting our numerical algorithm we describe how to numerically handle this problem of weak hyperbolicity, as well as how to choose an appropriate gauge condition. The resulting scheme is applied to several numerical test cases.

Scott A. Moe, James A. Rossmanith, David C. Seal

This work introduces a single-stage, single-step method for the compressible Euler equations that is provably positivity-preserving and can be applied on both Cartesian and unstructured meshes. This method is the first case of a single-stage, single-step method that is simultaneously high-order, positivity-preserving, and operates on unstructured meshes. Time-stepping is accomplished via the Lax-Wendroff approach, which is also sometimes called the Cauchy-Kovalevskaya procedure, where temporal derivatives in a Taylor series in time are exchanged for spatial derivatives. The Lax-Wendroff discontinuous Galerkin (LxW-DG) method developed in this work is formulated so that it looks like a forward Euler update but with a high-order time-extrapolated flux. In particular, the numerical flux used in this work is a linear combination of a low-order positivity-preserving contribution and a high-order component that can be damped to enforce positivity of the cell averages for the density and pressure for each time step. In addition to this flux limiter, a moment limiter is applied that forces positivity of the solution at finitely many quadrature points within each cell. The combination of the flux limiter and the moment limiter guarantees positivity of the cell averages from one time-step to the next. Finally, a simple shock capturing limiter that uses the same basic technology as the moment limiter is introduced in order to obtain non-oscillatory results. The resulting scheme can be extended to arbitrary order without increasing the size of the effective stencil. We present numerical results in one and two space dimensions that demonstrate the robustness of the proposed scheme.

Yongtao Cheng, James A. Rossmanith

Quadrature-based moment-closure methods are a class of approximations that replace high-dimensional kinetic descriptions with lower-dimensional fluid models. In this work we investigate some of the properties of a sub-class of these methods based on bi-delta, bi-Gaussian, and bi-B-spline representations. We develop a high-order discontinuous Galerkin (DG) scheme to solve the resulting fluid systems. Finally, via this high-order DG scheme and Strang operator splitting to handle the collision term, we simulate the fluid-closure models in the context of the Vlasov-Poisson-Fokker-Planck system in the high-field limit. We demonstrate numerically that the proposed scheme is asymptotic-preserving in the high-field limit.

Pierson T. Guthrey, James A. Rossmanith

Discontinuous Galerkin (DG) methods for hyperbolic partial differential equations (PDEs) with explicit time-stepping schemes, such as strong stability-preserving Runge-Kutta (SSP-RK), suffer from time-step restrictions that are significantly worse than what a simple Courant-Friedrichs-Lewy (CFL) argument requires. In particular, the maximum stable time-step scales inversely with the highest degree in the DG polynomial approximation space and becomes progressively smaller with each added spatial dimension. In this work we introduce a novel approach that we have dubbed the regionally implicit discontinuous Galerkin (RIDG) method to overcome these small time-step restrictions. The RIDG method is based on an extension of the Lax-Wendroff DG (LxW-DG) method, which previously had been shown to be equivalent to a predictor-corrector approach, where the predictor is a locally implicit spacetime method (i.e., the predictor is something like a block-Jacobi update for a fully implicit spacetime DG method). The corrector is an explicit method that uses the spacetime reconstructed solution from the predictor step. In this work we modify the predictor to include not just local information, but also neighboring information. With this modification we show that the stability is greatly enhanced; in particular, we show that we are able to remove the polynomial degree dependence of the maximum time-step and show how this extends to multiple spatial dimensions. A semi-analytic von Neumann analysis is presented to theoretically justify the stability claims. Convergence and efficiency studies for linear and nonlinear problems in multiple dimensions are accomplished using a MATLAB code that can be freely downloaded.

James A. Rossmanith

Residual distributions (RD) schemes are a class of of high-resolution finite volume methods for unstructured grids. A key feature of these schemes is that they make use of genuinely multidimensional (approximate) Riemann solvers as opposed to the piecemeal 1D Riemann solvers usually employed by finite volume methods. In 1D, LeVeque and Pelanti [J. Comp. Phys. 172, 572 (2001)] showed that many of the standard approximate Riemann solver methods (e.g., the Roe solver, HLL, Lax-Friedrichs) can be obtained from applying an exact Riemann solver to relaxation systems of the type introduced by Jin and Xin [Comm. Pure Appl. Math. 48, 235 (1995)]. In this work we extend LeVeque and Pelanti's results and obtain a multidimensional relaxation system from which multidimensional approximate Riemann solvers can be obtained. In particular, we show that with one choice of parameters the relaxation system yields the standard N-scheme. With another choice, the relaxation system yields a new Riemann solver, which can be viewed as a genuinely multidimensional extension of the local Lax-Friedrichs scheme. This new Riemann solver does not require the use Roe-Struijs-Deconinck averages, nor does it require the inversion of an m-by-m matrix in each computational grid cell, where $m$ is the number of conserved variables. Once this new scheme is established, we apply it on a few standard cases for the 2D compressible Euler equations of gas dynamics. We show that through the use of linear-preserving limiters, the new approach produces numerical solutions that are comparable in accuracy to the N-scheme, despite being computationally less expensive.

Christiane Helzel, James A. Rossmanith, Bertram Taetz

Numerical methods for solving the ideal magnetohydrodynamic (MHD) equations in more than one space dimension must confront the challenge of controlling errors in the discrete divergence of the magnetic field. One approach that has been shown successful in stabilizing MHD calculations are constrained transport (CT) schemes. CT schemes can be viewed as predictor-corrector methods for updating the magnetic field, where a magnetic field value is first predicted by a method that does not exactly preserve the divergence-free condition on the magnetic field, followed by a correction step that aims to control these divergence errors. In Helzel et al. (2011) the authors presented an unstaggered constrained transport method for the MHD equations on 3D Cartesian grids. In this work we generalize the method of Helzel et al. (2011) in three important ways: (1) we remove the need for operator splitting by switching to an appropriate method of lines discretization and coupling this with a non-conservative finite volume method for the magnetic vector potential equation, (2) we increase the spatial and temporal order of accuracy of the entire method to third order, and (3) we develop the method so that it is applicable on both Cartesian and logically rectangular mapped grids. The evolution equation for the magnetic vector potential is solved using a non-conservative finite volume method. The curl of the magnetic potential is computed via a third-order accurate discrete operator that is derived from appropriate application of the divergence theorem and subsequent numerical quadrature on element faces. Special artificial resistivity limiters are used to control unphysical oscillations in the magnetic potential and field components across shocks. Test computations are shown that confirm third order accuracy for smooth test problems and high-resolution for test problems with shock waves.

Evan Alexander Johnson, James A. Rossmanith

Numerical methods for hyperbolic conservation laws are needed that efficiently mimic the constraints satisfied by exact solutions, including material conservation and positivity, while also maintaining high-order accuracy and numerical stability. Discontinuous Galerkin (DG) and WENO schemes allow efficient high-order accuracy while maintaining conservation. Positivity limiters developed by Zhang and Shu ensure a minimum time step for which positivity of cell average quantities is maintained without sacrificing conservation or formal accuracy; this is achieved by linearly damping the deviation from the cell average just enough to enforce a cell positivity condition that requires positivity at boundary nodes and strategically chosen interior points. We assume that the set of positive states is convex; it follows that positivity is equivalent to scalar positivity of a collection of affine functionals. Based on this observation, we generalize the method of Zhang and Shu to a framework that we call outflow positivity limiting: First, enforce positivity at boundary nodes. If wave speed desingularization is needed, cap wave speeds at physically justified maxima by using remapped states to calculate fluxes. Second, apply linear damping again to cap the boundary average of all positivity functionals at the maximum possible (relative to the cell average) for a scalar-valued representation positive in each mesh cell. This be done by enforcing positivity of the retentional, an affine combination of the cell average and the boundary average, in the same way that Zhang and Shu would enforce positivity at a single point (and with similar computational expense). Third, limit the time step so that cell outflow is less than the initial cell content. This framework guarantees essentially the same positivity-preserving time step as is guaranteed if positivity is enforced at every point in the mesh cell.

E. Alec Johnson, James A. Rossmanith

For the first time to our knowledge, we demonstrate fast magnetic reconnection near a magnetic null point in a fluid model of collisionless pair plasma without resorting to the contrivance of anomalous resistivity. In particular, we demonstrate that fast reconnection occurs in an anisotropic adiabatic two-fluid model of collisionless pair plasma with relaxation toward isotropy for a broad range of isotropization rates. For very rapid isotropization we see fast reconnection, but instabilities eventually arise that cause numerical error and cast doubt on the simulated behavior.

Scott A. Moe, James A. Rossmanith, David C. Seal

The discontinuous Galerkin (DG) finite element method when applied to hyperbolic conservation laws requires the use of shock-capturing limiters in order to suppress unphysical oscillations near large solution gradients. In this work we develop a novel shock-capturing limiter that combines key ideas from the limiter of Barth and Jespersen [AIAA-89-0366 (1989)] and the maximum principle preserving (MPP) framework of Zhang and Shu [Proc. R. Soc. A, 467 (2011), pp. 2752--2776]. The limiting strategy is based on traversing the mesh element-by-element in order to (1) find local upper and lower bounds on user-defined variables by sampling these variables on neighboring elements, and (2) to then enforce these local bounds by minimally damping the high-order corrections. The main advantages of this limiting strategy is that it is simple to implement, effective at shock capturing, and retains high-order accuracy of the solution in smooth regimes. The resulting numerical scheme is applied to several standard numerical tests in both one and two-dimensions and on both Cartesian and unstructured grids. These tests are used as benchmarks to verify and assess the accuracy and robustness of the method.