A Maximum-Principle-Satisfying High-Order Finite Volume Compact WENO Scheme for Scalar Conservation Laws with Applications in Incompressible Flows

In this paper, a maximum-principle-satisfying finite volume compact scheme is proposed for solving scalar hyperbolic conservation laws. The scheme combines weighted essentially non-oscillatory schemes (WENO) with a class of compact schemes under a finite volume framework, in which the nonlinear WENO weights are coupled with lower order compact stencils. The maximum-principle-satisfying polynomial rescaling limiter in Zhang and Shu (J Comput Phys 229:3091–3120, 2010, Proc R Soc A Math Phys Eng Sci 467:2752–2776, 2011) is adopted to construct the present schemes at each stage of an explicit Runge–Kutta method, without destroying high order accuracy and conservativity. Numerical examples for one and two dimensional problems including incompressible flows are presented to assess the good performance, maximum principle preserving, essentially non-oscillatory and high resolution of the proposed method.     

Finite Volume HWENO Schemes for Nonconvex Conservation Laws

Following the previous work of Qiu and Shu (SIAM J Sci Comput 31: 584–607, 2008), we investigate the performance of Hermite weighted essentially non-oscillatory (HWENO) scheme for nonconvex conservation laws. Similar to many other high order methods, we show that the finite volume HWENO scheme performs poorly for some nonconvex conservation laws. We modify the scheme around the nonconvex regions, based on a first order monotone scheme and a second entropic projection, to ensure entropic convergence. Extensive numerical tests are performed. Compare with the earlier work of Qiu and Shu which focuses on 1D scalar problems, we apply the modified schemes (both WENO and HWENO) to one-dimensional Euler system with nonconvex equation of state and two-dimensional problems.     

Superconvergent Two-Grid Methods for Elliptic Eigenvalue Problems

Some numerical algorithms for elliptic eigenvalue problems are proposed, analyzed, and numerically tested. The methods combine advantages of the two-grid algorithm (Xu and Zhou in Math Comput 70(233):17–25, 2001), the two-space method (Racheva and Andreev in Comput Methods Appl Math 2:171–185, 2002), the shifted inverse power method (Hu and Cheng in Math Comput 80:1287–1301, 2011; Yang and Bi in SIAM J Numer Anal 49:1602–1624, 2011), and the polynomial preserving recovery enhancing technique (Naga et al. in SIAM J Sci Comput 28:1289–1300, 2006). Our new algorithms compare favorably with some existing methods and enjoy superconvergence property.     

Direct Discontinuous Galerkin Method and Its Variations for Second Order Elliptic Equations

In this paper, we study direct discontinuous Galerkin method (Liu and Yan in SIAM J Numer Anal 47(1):475–698, 2009) and its variations (Liu and Yan in Commun Comput Phys 8(3):541–564, 2010; Vidden and Yan in J Comput Math 31(6):638–662, 2013; Yan in J Sci Comput 54(2–3):663–683, 2013) for 2nd order elliptic problems. A priori error estimate under energy norm is established for all four methods. Optimal error estimate under \(L^2\) norm is obtained for DDG method with interface correction (Liu and Yan in Commun Comput Phys 8(3):541–564, 2010) and symmetric DDG method (Vidden and Yan in J Comput Math 31(6):638–662, 2013). A series of numerical examples are carried out to illustrate the accuracy and capability of the schemes. Numerically we obtain optimal \((k+1)\)th order convergence for DDG method with interface correction and symmetric DDG method on nonuniform and unstructured triangular meshes. An interface problem with discontinuous diffusion coefficients is investigated and optimal \((k+1)\)th order accuracy is obtained. Peak solutions with sharp transitions are captured well. Highly oscillatory wave solutions of Helmholz equation are well resolved.     

Hybrid Well-balanced WENO Schemes with Different Indicators for Shallow Water Equations

In (J. Comput. Phys. 229: 8105–8129, 2010), Li and Qiu investigated the hybrid weighted essentially non-oscillatory (WENO) schemes with different indicators for Euler equations of gas dynamics. In this continuation paper, we extend the method to solve the one- and two-dimensional shallow water equations with source term due to the non-flat bottom topography, with a goal of obtaining the same advantages of the schemes for the Euler equations, such as the saving computational cost, essentially non-oscillatory property for general solution with discontinuities, and the sharp shock transition. Extensive simulations in one- and two-dimensions are provided to illustrate the behavior of this procedure.     

High Order Positivity- and Bound-Preserving Hybrid Compact-WENO Finite Difference Scheme for the Compressible Euler Equations

Based on the same hybridization framework of Don et al. (SIAM J Sci Comput 38:A691–A711 2016), an improved hybrid scheme employing the nonlinear 5th-order characteristic-wise WENO-Z5 finite difference scheme for capturing high gradients and discontinuities in an essentially non-oscillatory manner and the linear 5th-order conservative compact upwind (CUW5) scheme for resolving the fine scale structures in the smooth regions of the solution in an efficient and accurate manner is developed. By replacing the 6th-order non-dissipative compact central scheme (CCD6) with the CUW5 scheme, which has a build-in dissipation, there is no need to employ an extra high order smoothing procedure to mitigate any numerical oscillations that might appear in an hybrid scheme. The high order multi-resolution algorithm of Harten is employed to detect the smoothness of the solution. To handle the problems with extreme conditions, such as high pressure and density ratios and near vacuum states, and detonation diffraction problems, we design a positivity- and bound-preserving limiter by extending the one developed in Hu et al. (J Comput Phys 242, 2013) for solving the high Mach number jet flows, detonation diffraction problems and detonation passing multiple obstacles problems. Extensive one- and two-dimensional shocked flow problems demonstrate that the new hybrid scheme is less dispersive and less dissipative, and allows a potential speedup up to a factor of more than one and half times faster than the WENO-Z5 scheme.     

WSGD-OSC Scheme for Two-Dimensional Distributed Order Fractional Reaction–Diffusion Equation

In this paper, a new numerical approximation is discussed for the two-dimensional distributed-order time fractional reaction–diffusion equation. Combining with the idea of weighted and shifted Grünwald difference (WSGD) approximation (Tian et al. in Math Comput 84:1703–1727, 2015; Wang and Vong in J Comput Phys 277:1–15, 2014) in time, we establish orthogonal spline collocation (OSC) method in space. A detailed analysis shows that the proposed scheme is unconditionally stable and convergent with the convergence order \(\mathscr {O}(\tau ^2+\Delta \alpha ^2+h^{r+1})\), where \(\tau , \Delta \alpha , h\) and r are, respectively the time step size, step size in distributed-order variable, space step size, and polynomial degree of space. Interestingly, we prove that the proposed WSGD-OSC scheme converges with the second-order in time, where OSC schemes proposed previously (Fairweather et al. in J Sci Comput 65:1217–1239, 2015; Yang et al. in J Comput Phys 256:824–837, 2014) can at most achieve temporal accuracy of order which depends on the order of fractional derivatives in the equations and is usually less than two. Some numerical results are also given to confirm our theoretical prediction.     

Weights Design For Maximal Order WENO Schemes

Weighted essentially non-oscillatory (WENO) finite difference schemes, developed by Liu et al. (Comput Phys 115(1):200–212, 1994) and improved by Jiang and Shu (Comput Phys 126(1):202–228, 1996), are one of the most popular methods to approximate the solutions of hyperbolic equations. But these schemes fail to provide maximal order accuracy near smooth extrema, where the first derivative of the solution becomes zero. Some authors have addressed this problem with different weight designs. In this paper we focus on the weights proposed by Yamaleev and Carpenter (J Comput Phys 228:4248–4272, 2009). They propose new weights to provide faster weight convergence than those presented in Borges et al. (J Comput Phys 227:3191–3211, 2008) and deduce some constraints on the weights parameters to guarantee that the WENO scheme has maximal order for sufficiently smooth solutions with an arbitrary number of vanishing derivatives. We analyze the scheme with the weights proposed in Yamaleev and Carpenter (J Comput Phys 228:4248–4272, 2009) and prove that near discontinuities it achieves worse orders than classical WENO schemes. In order to solve these accuracy problems, we define new weights, based on those proposed in Yamaleev and Carpenter (J Comput Phys 228:4248–4272, 2009), and get some constraints on the weights parameters to guarantee maximal order accuracy for the resulting schemes.     

Flux-based level-set method for two-phase flows on unstructured grids

In this paper we deal with the application of the flux-based level set method to moving interface computations on unstructured grids. The focus lies on the overcoming of the known difficulties of level set methods, e.g. accurate computations of important geometric properties, reliable and precise reinitialization of the level set function and the adaption of standard discretization methods to the moving boundary case. The basic building block of our approach is the high-resolution flux-based level set method for general advection equation (Frolkovi? and Mikula in SIAM J Sci Comput 29(2):579–597, 2007, Frolkovi? and Wehner in Comput Vis Sci 12(6):626–650, 2009). We extend this method for the problem of reinitialization of the level set function on unstructured grids by using quadratic interpolation to compute distances for nodes close to the interface. To realize numerical simulation for some applications with moving boundaries, we adapt the approach of ghost fluid method (Gibou and Fedkiw in J Comput Phys 202:577–601, 2005) for unstructured grids. The idea is to describe the development of the moving boundary with a level set formulation while the computational grid remains fixed and the boundary conditions are enforced using some extrapolation. Our main motivation is the numerical solution of two-phase incompressible flow problems. Additionally to previously mentioned steps, we introduce further numerical schemes in the framework of finite volume discretization for the flow. Possible jumps of the pressure and the directional derivative of velocity at the interface are modeled directly within the method using the approach of extended approximation spaces. Besides that, an algorithm for the computations of curvature is considered that exhibits the second order accuracy for some examples. Numerical experiments are provided for the presented methods.     

Modeling and Simulation of Laser Processing of Particulate-Functionalized Materials

The objective of this paper is to focus on one of the “building blocks” of additive manufacturing technologies, namely selective laser-processing of particle-functionalized materials. Following a series of work in Zohdi (Int J Numer Methods Eng 53:1511–1532, 2002; Philos Trans R Soc Math Phys Eng Sci 361(1806):1021–1043, 2003; Comput Methods Appl Mech Eng 193(6–8):679–699, 2004; Comput Methods Appl Mech Eng 196:3927–3950, 2007; Int J Numer Methods Eng 76:1250–1279, 2008; Comput Methods Appl Mech Eng 199:79–101, 2010; Arch Comput Methods Eng 1–17. doi: 10.1007/s11831-013-9092-6, 2013; Comput Mech Eng Sci 98(3):261–277, 2014; Comput Mech 54:171–191, 2014; J Manuf Sci Eng ASME doi: 10.1115/1.4029327, 2015; CIRP J Manuf Sci Technol 10:77–83, 2015; Comput Mech 56:613–630, 2015; Introduction to computational micromechanics. Springer, Berlin, 2008; Introduction to the modeling and simulation of particulate flows. SIAM (Society for Industrial and Applied Mathematics), Philadelphia, 2007; Electromagnetic properties of multiphase dielectrics: a primer on modeling, theory and computation. Springer, Berlin, 2012), a laser-penetration model, in conjunction with a Finite Difference Time Domain Method using an immersed microstructure method, is developed. Because optical, thermal and mechanical multifield coupling is present, a recursive, staggered, temporally-adaptive scheme is developed to resolve the internal microstructural fields. The time step adaptation allows the numerical scheme to iteratively resolve the changing physical fields by refining the time-steps during phases of the process when the system is undergoing large changes on a relatively small time-scale and can also enlarge the time-steps when the processes are relatively slow. The spatial discretization grids are uniform and dense enough to capture fine-scale changes in the fields. The microstructure is embedded into the spatial discretization and the regular grid allows one to generate a matrix-free iterative formulation which is amenable to rapid computation, with minimal memory requirements, making it ideal for laptop computation. Numerical examples are provided to illustrate the modeling and simulation approach, which by design, is straightforward to computationally implement, in order to be easily utilized by researchers in the field. More advanced conduction models, based on thermal-relaxation, which are a key feature of fast-pulsing laser technologies, are also discussed.     

Optimal Spectral Schemes Based on Generalized Prolate Spheroidal Wave Functions of Order $$$$

We introduce a family of generalized prolate spheroidal wave functions (PSWFs) of order \(-1,\) and develop new spectral schemes for second-order boundary value problems. Our technique differs from the differentiation approach based on PSWFs of order zero in Kong and Rokhlin (Appl Comput Harmon Anal 33(2):226–260, 2012); in particular, our orthogonal basis can naturally include homogeneous boundary conditions without the re-orthogonalization of Kong and Rokhlin (2012). More notably, it leads to diagonal systems or direct “explicit” solutions to 1D Helmholtz problems in various situations. Using a rule optimally pairing the bandwidth parameter and the number of basis functions as in Kong and Rokhlin (2012), we demonstrate that the new method significantly outperforms the Legendre spectral method in approximating highly oscillatory solutions. We also conduct a rigorous error analysis of this new scheme. The idea and analysis can be extended to generalized PSWFs of negative integer order for higher-order boundary value and eigenvalue problems.     

Improvement of Convergence to Steady State Solutions of Euler Equations with the WENO Schemes

The convergence to steady state solutions of the Euler equations for high order weighted essentially non-oscillatory (WENO) finite difference schemes with the Lax-Friedrichs flux splitting (Jiang and Shu, in J. Comput. Phys. 126:202–228, 1996) is investigated. Numerical evidence in Zhang and Shu (J. Sci. Comput. 31:273–305, 2007) indicates that there exist slight post-shock oscillations when we use high order WENO schemes to solve problems containing shock waves. Even though these oscillations are small in their magnitude and do not affect the “essentially non-oscillatory” property of the WENO schemes, they are indeed responsible for the numerical residue to hang at the truncation error level of the scheme instead of settling down to machine zero. Differently from the strategy adopted in Zhang and Shu (J. Sci. Comput. 31:273–305, 2007), in which a new smoothness indicator was introduced to facilitate convergence to steady states, in this paper we study the effect of the local characteristic decomposition on steady state convergence. Numerical tests indicate that the slight post-shock oscillation has a close relationship with the local characteristic decomposition process. When this process is based on an average Jacobian at the cell interface using the Roe average, as is the standard procedure for WENO schemes, such post-shock oscillation appears. If we instead use upwind-biased interpolation to approximate the physical variables including the velocity and enthalpy on the cell interface to compute the left and right eigenvectors of the Jacobian for the local characteristic decomposition, the slight post-shock oscillation can be removed or reduced significantly and the numerical residue settles down to lower values than other WENO schemes and can reach machine zero for many test cases. This new procedure is also effective for higher order WENO schemes and for WENO schemes with different smoothness indicators.     

A modified fifth-order WENO scheme for hyperbolic conservation laws

Recently a WENO scheme, with smoothness indicators constructed based on $L_1$ measure is introduced by Ha et al. (2013) and the improved version of this scheme is presented by Kim et al. (2016), referred to as WENO-NS and WENO-P schemes respectively. These schemes perform better than the existing many fifth-order WENO schemes for the problems which contain discontinuities and attain fifth-order accuracy at the critical points where the first derivative vanishes but not at the points where the second derivatives are zero. This paper deals with modification of the above said methods to obtain a new fifth-order weighted essentially non-oscillatory (WENO) scheme. A new global-smoothness indicator is proposed which shows an improved behavior over the solutions of WENO-NS and WENO-P schemes and the proposed scheme attains an optimal order of approximation, even at the critical points where the first and second derivatives vanish but not the third derivative. Examples are taken in the numeric section to check the robustness and accuracy of the proposed scheme for one and two-dimensional system of Euler equations.     

Vertex-Centered Linearity-Preserving Schemes for Nonlinear Parabolic Problems on Polygonal Grids

On arbitrary polygonal grids, a family of vertex-centered finite volume schemes are suggested for the numerical solution of the strongly nonlinear parabolic equations arising in radiation hydrodynamics and magnetohydrodynamics. We define the primary unknowns at the cell vertices and derive the schemes along the linearity-preserving approach. Since we adopt the same cell-centered diffusion coefficients as those in most existing finite volume schemes, it is required to introduce some auxiliary unknowns at the cell centers in the case of nonlinear diffusion coefficients. A second-order positivity-preserving algorithm is then suggested to interpolate these auxiliary unknowns via the primary ones. All the schemes lead to symmetric and positive definite linear systems and their stability can be rigorously analyzed under some standard and weak geometry assumptions. More interesting is that these vertex-centered schemes do not have the so-called numerical heat-barrier issue suffered by many existing cell-centered or hybrid schemes (Lipnikov et al. in J Comput Phys 305:111–126, 2016). Numerical experiments are also presented to show the efficiency and robustness of the schemes in simulating nonlinear parabolic problems.     

New WENO Smoothness Indicators Computationally Efficient in the Presence of Corner Discontinuities

This paper is devoted to the construction and analysis of a new smoothness index for WENO interpolation capable of dealing with corner discontinuities. The new smoothness index presented is initially developed for the point-value framework of Harten’s multiresolution. Even so, the ideas about how to extend the results to the cell-average framework are presented. The new smoothness index is inspired by the one proposed in Jiang and Shu (J Comput Phys 126(1):202–228, 1996). This index works very well for jump discontinuities as it was originally designed for the context of conservation laws in order to deal with problems that contain shocks and complicated fluid-structure interactions. Even so, it is easy to check that the mentioned index does not provide an appropriate performance for corner discontinuities. Our aim is to rise the order of accuracy of WENO interpolation near corner discontinuities. In order to do so, we will modify the original smoothness index proposed by Jiang and Shu such that the discontinuities in the first derivative of the function contribute effectively to the index. The modification proposed will produce a variation in the weights of WENO when dealing with a corner, that do not appear when using the smoothness indexes proposed by Jiang and Shu. The variation in the weights induced by the modification of the smoothness index will allow adaption to corner discontinuities, maintaining the adaption to jumps provided by the original smoothness index proposed by Jiang and Shu. The strategy proposed in Aràndiga et al. (SIAM J Numer Anal 49(2):893–915, 2011) can be adapted such that the accuracy is maintained near critical points at smooth zones.     

High Order Multi-dimensional Characteristics Tracing for the Incompressible Euler Equation and the Guiding-Center Vlasov Equation

In this paper, we propose a high order characteristics tracing scheme for the two-dimensional nonlinear incompressible Euler system in vorticity stream function formulation and the guiding center Vlasov model. Such a scheme is incorporated into a semi-Lagrangian finite difference WENO framework for simulating the aforementioned model equations. This is an extension of our earlier work on high order characteristics tracing scheme for the 1D nonlinear Vlasov–Poisson system (Qiu and Russo in J Sci Comput 71:414–434, 2017). The effectiveness of the proposed scheme is demonstrated numerically by an extensive set of test cases.     

A Multi-Level Mixed Element Method for the Eigenvalue Problem of Biharmonic Equation

In this paper, we discuss approximating the eigenvalue problem of biharmonic equation. We first present an equivalent mixed formulation which admits natural nested discretization. Then, we present multi-level finite element schemes by implementing the algorithm as in Lin and Xie (Math Comput 84:71–88, 2015) to the nested discretizations on a series of nested grids. The multi-level mixed scheme for the biharmonic eigenvalue problem possesses optimal convergence rate and optimal computational cost. Both theoretical analysis and numerical verifications are presented.     

Mesh generation for confined fusion plasma simulation

XGC1 and M3D-C ¹ are two fusion plasma simulation codes being developed at Princeton Plasma Physics Laboratory. XGC1 uses the particle-in-cell method to simulate gyrokinetic neoclassical physics and turbulence (Chang et al. Phys Plasmas 16(5):056108, 2009; Ku et al. Nucl Fusion 49:115021, 2009; Admas et al. J Phys 180(1):012036, 2009). M3D-\(C^1\) solves the two-fluid resistive magnetohydrodynamic equations with the \(C^1\) finite elements (Jardin J comput phys 200(1):133–152, 2004; Jardin et al. J comput Phys 226(2):2146–2174, 2007; Ferraro and Jardin J comput Phys 228(20):7742–7770, 2009; Jardin J comput Phys 231(3):832–838, 2012; Jardin et al. Comput Sci Discov 5(1):014002, 2012; Ferraro et al. Sci Discov Adv Comput, 2012; Ferraro et al. International sherwood fusion theory conference, 2014). This paper presents the software tools and libraries that were combined to form the geometry and automatic meshing procedures for these codes. Specific consideration has been given to satisfy the mesh configuration and element shape quality constraints of XGC1 and M3D-\(C^1\).