Search

The quadrature-based method of moments (QMOM) offers a promising class of approximation techniques for reducing kinetic equations to fluid equations that are valid beyond thermodynamic equilibrium. In this work, we study a particular five-moment variant of QMOM known as HyQMOM and establish that this system is moment-invertible over a convex region in solution space. We then develop a high-order discontinuous Galerkin (DG) scheme for solving the resulting fluid system. The scheme is based on a predictor-corrector approach, where the prediction is a localized space-time DG scheme. The nonlinear algebraic system in this prediction is solved using a Picard iteration. The correction is a straightforward explicit update based on the time-integral of the evolution equation, where the space-time prediction replaces all instances of the exact solution. In the absence of limiters, the high-order scheme does not guarantee that solutions remain in the convex set over which HyQMOM is moment-realizable. To overcome this, we introduce novel limiters that rigorously guarantee that the computed solution does not leave the convex set of realizable solutions, thus guaranteeing the hyperbolicity of the system. We develop positivity-preserving limiters in both the prediction and correction steps and an oscillation limiter that damps unphysical oscillations near shocks. We also develop a novel extension of this scheme to include a BGK collision operator; the proposed method is shown to be asymptotic-preserving in the high-collision limit. The HyQMOM and the HyQMOM-BGK solvers are verified on several test cases, demonstrating high-order accuracy on smooth problems and shock-capturing capability on problems with shocks. The asymptotic-preserving property of the HyQMOM-BGK solver is also numerically verified., Comment: 54 pages, 9 figures, 4 tables

The Cahn-Hilliard Navier-Stokes (CHNS) system provides a computationally tractable model that can be used to effectively capture interfacial dynamics in two-phase fluid flows. In this work, we present a semi-implicit, projection-based finite element framework for solving the CHNS system. We use a projection-based semi-implicit time discretization for the Navier-Stokes equation and a fully-implicit time discretization for the Cahn-Hilliard equation. We use a conforming continuous Galerkin (cG) finite element method in space equipped with a residual-based variational multiscale (RBVMS) formulation. Pressure is decoupled using a projection step, which results in two linear positive semi-definite systems for velocity and pressure, instead of the saddle point system of a pressure-stabilized method. All the linear systems are solved using an efficient and scalable algebraic multigrid (AMG) method. We deploy this approach on a massively parallel numerical implementation using parallel octree-based adaptive meshes. The overall approach allows the use of relatively large time steps with much faster time-to-solve than similar fully-implicit methods. We present comprehensive numerical experiments showing detailed comparisons with results from the literature for canonical cases, including the single bubble rise and Rayleigh-Taylor instability., Comment: 51 pages, 18 Figures, version accepted for publication Journal of Computational Physics. arXiv admin note: text overlap with arXiv:2009.06628

The thermal radiative transfer (TRT) equations form an integro-differential system that describes the propagation and collisional interactions of photons. Computing accurate and efficient numerical solutions TRT are challenging for several reasons, the first of which is that TRT is defined on a high-dimensional phase. In order to reduce the dimensionality of the phase space, classical approaches such as the P$_N$ (spherical harmonics) or the S$_N$ (discrete ordinates) ansatz are often used in the literature. In this work, we introduce a novel approach: the hybrid discrete (H$^T_N$) approximation to the radiative thermal transfer equations. This approach acquires desirable properties of both P$_N$ and S$_N$, and indeed reduces to each of these approximations in various limits: H$^1_N$ $\equiv$ P$_N$ and H$^T_0$ $\equiv$ S$_T$. We prove that H$^T_N$ results in a system of hyperbolic equations for all $T\ge 1$ and $N\ge 0$. Another challenge in solving the TRT system is the inherent stiffness due to the large timescale separation between propagation and collisions, especially in the diffusive (i.e., highly collisional) regime. This stiffness challenge can be partially overcome via implicit time integration, although fully implicit methods may become computationally expensive due to the strong nonlinearity and system size. On the other hand, explicit time-stepping schemes that are not also asymptotic-preserving in the highly collisional limit require resolving the mean-free path between collisions, making such schemes prohibitively expensive. In this work we develop a numerical method that is based on a nodal discontinuous Galerkin discretization in space, coupled with a semi-implicit discretization in time. We conduct several numerical experiments to verify the accuracy, efficiency, and robustness of the H$^T_N$ ansatz and the numerical discretizations., Comment: 33 pages, 11 figures, 1 table

In this work we investigate the parallel scalability of the numerical method developed in Guthrey and Rossmanith [The regionally implicit discontinuous Galerkin method: Improving the stability of DG-FEM, SIAM J. Numer. Anal. (2019)]. We develop an implementation of the regionally-implicit discontinuous Galerkin (RIDG) method in DoGPack, which is an open source C++ software package for discontinuous Galerkin methods. Specifically, we develop and test a hybrid OpenMP and MPI parallelized implementation of DoGPack with the goal of exploring the efficiency and scalability of RIDG in comparison to the popular strong stability-preserving Runge-Kutta discontinuous Galerkin (SSP-RKDG) method. We demonstrate that RIDG methods are able to hide communication latency associated with distributed memory parallelism, due to the fact that almost all of the work involved in the method is highly localized to each element, producing a localized prediction for each region. We demonstrate the enhanced efficiency and scalability of the of the RIDG method and compare it to SSP-RKDG methods and show extensibility to very high order schemes. The two-dimensional scaling study is performed on machines at the Institute for Cyber-Enabled Research at Michigan State University, using up to 1440 total cores on Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz CPUs. The three dimensional scaling study is performed on Livermore Computing clusters at at Lawrence Livermore National Laboratory, using up to 28672 total cores on Intel Xeon CLX-8276L CPUs with Omni-Path interconnects., Comment: 26 pages, 2 figures, 6 tables

We present a fully-coupled, implicit-in-time framework for solving a thermodynamically-consistent Cahn-Hilliard Navier-Stokes system that models two-phase flows. In this work, we extend the block iterative method presented in Khanwale et al. [\textit{Simulating two-phase flows with thermodynamically consistent energy stable Cahn-Hilliard Navier-Stokes equations on parallel adaptive octree based meshes}, J. Comput. Phys. (2020)], to a fully-coupled, provably second-order accurate scheme in time, while maintaining energy-stability. The new method requires fewer matrix assemblies in each Newton iteration resulting in faster solution time. The method is based on a fully-implicit Crank-Nicolson scheme in time and a pressure stabilization for an equal order Galerkin formulation. That is, we use a conforming continuous Galerkin (cG) finite element method in space equipped with a residual-based variational multiscale (RBVMS) procedure to stabilize the pressure. We deploy this approach on a massively parallel numerical implementation using parallel octree-based adaptive meshes. We present comprehensive numerical experiments showing detailed comparisons with results from the literature for canonical cases, including the single bubble rise, Rayleigh-Taylor instability, and lid-driven cavity flow problems. We analyze in detail the scaling of our numerical implementation., Comment: 55 pages, 18 figures, 3 tables, version accepted in Computer Physics Communications

We report on simulations of two-phase flows with deforming interfaces at various density contrasts by solving thermodynamically consistent Cahn-Hilliard Navier-Stokes equations. An (essentially) unconditionally energy-stable Crank-Nicolson-type time integration scheme is used. Detailed proofs of energy stability of the semi-discrete scheme and for the existence of solutions of the advective-diffusive Cahn-Hilliard operator are provided. In space we discretize with a conforming continuous Galerkin finite element method in conjunction with a residual-based variational multi-scale (VMS) approach in order to provide pressure stabilization. We deploy this approach on a massively parallel numerical implementation using fast octree-based adaptive meshes. A detailed scaling analysis of the solver is presented. Numerical experiments showing convergence and validation with experimental results from the literature are presented for a large range of density ratios., Comment: 59 pages, 22 figures, submitted to Journal of Computational Physics; Version accepted by the journal

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., Comment: 37 pages, 7 figures, 4 tables, 5 algorithms

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., Comment: 26 pages, 4 figures, 8 tables

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., Comment: 28 pages, 9 figures

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., Comment: 20 pages

The modification of the celebrated Yee scheme from Maxwell equations to magnetohydrodynamics is often referred to as the constrained transport approach. Constrained transport can be viewed as a sort of predictor-corrector method for updating the magnetic field, where a magnetic field value is first predicted by a method that does not preserve the divergence-free condition on the magnetic field, followed by a correction step that aims to control these divergence errors. This strategy has been successfully used in conjunction with a variety of shock-capturing methods including WENO, central, and wave propagation schemes. In this work we show how to extend the basic CT framework to the discontinuous Galerkin finite element method on both 2D and 3D Cartesian grids. We first review the entropy-stability theory for semi-discrete DG discretizations of ideal MHD, which rigorously establishes the need for a magnetic field that satisfies the following conditions: (1) the divergence of the magnetic field is zero on each element, and (2) the normal components of the magnetic field are continuous across element edges/faces. In order to achieve such a globally divergence-free magnetic field, we introduce a novel CT scheme that is based on two ingredients: (1) we introduce an element-centered magnetic vector potential that is updated via a discontinuous Galerkin scheme on the induction equation; and (2) we define a mapping that takes element-centered magnetic field values and element-centered magnetic vector potential values and creates on each edge /face a representation of the normal component of the magnetic field; this representation is then mapped back to the elements to create a globally divergence-free element-centered representation of the magnetic field. For problems with shock waves, we make use of so-called moment-based limiters to control oscillations in the conserved quantities., Comment: 26 pages, 6 figures

In this work we develop a class of high-order finite difference weighted essentially non-oscillatory (FD-WENO) schemes for solving the ideal magnetohydrodynamic (MHD) equations in 2D and 3D. The philosophy of this work is to use efficient high-order WENO spatial discretizations with high-order strong stability-preserving Runge-Kutta (SSP-RK) time-stepping schemes. Numerical results have shown that with such methods we are able to resolve solution structures that are only visible at much higher grid resolutions with lower-order schemes. The key challenge in applying such methods to ideal MHD is to control divergence errors in the magnetic field. We achieve this by augmenting the base scheme with a novel high-order constrained transport approach that updates the magnetic vector potential. The predicted magnetic field from the base scheme is replaced by a divergence-free magnetic field that is obtained from the curl of this magnetic potential. The non-conservative weakly hyperbolic system that the magnetic vector potential satisfies is solved using a version of FD-WENO developed for Hamilton-Jacobi equations. The resulting numerical method is endowed with several important properties: (1) all quantities, including all components of the magnetic field and magnetic potential, are treated as point values on the same mesh (i.e., there is no mesh staggering); (2) both the spatial and temporal orders of accuracy are fourth-order; (3) no spatial integration or multidimensional reconstructions are needed in any step; and (4) special limiters in the magnetic vector potential update are used to control unphysical oscillations in the magnetic field. Several 2D and 3D numerical examples are presented to verify the order of accuracy on smooth test problems and to show high-resolution on test problems that involve shocks., Comment: 39 pages, 9 figures, 4 tables

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., Comment: 32 pages, 15 figures

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., Comment: 24 pages, 5 figures

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., Comment: 29 pages, 5 figures, 3 tables

We simulate magnetic reconnection in the GEM problem using a two-fluid model with 10 moments for the electron fluid as well as the proton fluid. We show that use of 10 moments for electrons gives good qualitative agreement with the the electron pressure tensor components in published kinetic simulations., Comment: 13th International Conference on Hyperbolic Problems, June 2010

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.

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., Comment: 46 pages, 12 figures

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., Comment: 46 pages, 14 figures