Lax–Friedrichs Multigrid Fast Sweeping Methods for Steady State Problems for Hyperbolic Conservation Laws

A fast marching algorithm for the factored eikonal equation

In [F. Li, C.-W. Shu, Y.-T. Zhang, H. Zhao, J. Comput. Phys., 227 (2008) pp. 8191–8208], we developed a fast sweeping method based on a hybrid local solver which is a combination of a discontinuous Galerkin (DG) finite element solver and a first order finite difference solver for Eikonal equations. The method has second order accuracy in the $L^1$ norm and a very fast convergence speed, but only first order accuracy in the $L^\infty$ norm for the general cases. This is an obstacle to the design of higher order DG fast sweeping methods. In this paper, we overcome this problem by developing uniformly accurate DG fast sweeping methods for solving Eikonal equations. We design novel causality indicators which guide the information flow directions for the DG local solver. The values of these indicators are initially provided by the first order finite difference fast sweeping method, and they are updated during iterations along with the solution. We observe both a uniform second order accuracy in the $L^\infty$ norm (in smooth regions) and the fast convergence speed (linear computational complexity) in the numerical examples.

Hybrid massively parallel fast sweeping method for static Hamilton–Jacobi equations

Efficient and accurate traveltime calculations of seismic waves have important applications in traveltime tomography for initial velocity model building, Kirchhoff depth migration, earthquake location, etc. Because anisotropy significantly affects the traveltimes of seismic waves, neglecting it would result in inaccurate imaging and inversion of the underground structure. Therefore, it is necessary to consider anisotropy in the traveltime calculations. The fast sweeping method (FSM) has important applications in computing the anisotropic first-arrival traveltime. The conventional method, which solves a transformed traveltime quartic equation combined with the FSM, is well suited for general anisotropic media and does not rely on approximations to Christoffel’s equation. However, this method has the following problems: (1) the form of the traveltime quartic equation is highly complex and exhibits multiple solutions, requiring additional effort to find solution intervals and determine which solution meets the criteria, and (2) Solving high-order polynomial equations may encounter numerical instability, especially when coefficients undergo minor variations or rounding errors occur for 3D problems. To address the stability and efficiency issues in anisotropic traveltime calculation, we analyze the characteristics of the constructed triangular-pyramid local solver and the quartic coupled slowness equation for the qP and qSV waves in tilted transversely isotropic (TTI) media. We observe that the decomposed slowness equation yields only one or no solution in the triangular-pyramid local solver and satisfies monotonicity conditions. Therefore, we develop the use of the Newton method to solve the factorized slowness equation efficiently, thus addressing stability and computational efficiency concerns inherent in the conventional approach. For the qSH wave, its slowness equation is quadratic and simple to solve. Our method provides an efficient and stable procedure for the traveltime calculations of qP, qSV, and qSH waves in 3D general TTI media. Numerical examples verify the efficiency and accuracy of our method.

This paper presents a fast sweeping method (FSM) to calculate the first‐arrival traveltimes of the qP, qSV, and qSH waves in two‐dimensional (2D) transversely isotropic media, whose symmetry axis may have an arbitrary orientation (tilted transverse isotropy [TTI]). The method discretizes the anisotropic eikonal equation with finite difference approximations on a rectangular mesh and solves the discretized system iteratively with the Gauss‐Seidel iterations along alternating sweeping orderings. At each mesh point, a highly nonlinear equation is solved to update the numerical solution until its convergence. For solving the nonlinear equation, an interval that contains the solutions is first determined and partitioned into few subintervals such that each subinterval contains one solution; then, the false position method is applied on these subintervals to compute the solutions; after that, among all possible solutions for the discretized equation, a causality condition is imposed, and the minimum solution satisfying the causality condition is chosen to update the solution. For problems with a point‐source condition, the FSM is extended for solving the anisotropic eikonal equation after a factorization technique is applied to resolve the source singularities, which yields clean first‐order accuracy. When dealing with the triplication of the qSV wave, solutions corresponding to the minimal group velocity are chosen such that continuous solutions are computed. The accuracy, efficiency, and capability of the proposed method are demonstrated with numerical experiments.

We further investigate the high order positivity-preserving discontinuous Galerkin (DG) methods for linear hyperbolic and radiative transfer equations developed in Yuan et al. (SIAM J Sci Comput 38:A2987–A3019, 2016). The DG methods in Yuan et al. (2016) can maintain positivity and high order accuracy, but they rely both on the scaling limiter in Zhang and Shu (J Comput Phys 229:8918–8934, 2010) and a rotational limiter, the latter may alter cell averages of the unmodulated DG scheme, thereby affecting conservation. Even though a Lax–Wendroff type theorem is proved in Yuan et al. (2016), guaranteeing convergence to weak solutions with correct shock speed when such rotational limiter is applied, it would still be desirable if a conservative DG method without changing the cell averages can be obtained which has both high order accuracy and positivity-preserving capability. In this paper, we develop and analyze such a DG method for both linear hyperbolic equations and radiative transfer equations. In the one-dimensional case, the method uses traditional DG space $$P^k$$ of piecewise polynomials of degree at most k. A key result is proved that the unmodulated DG solver in this case can maintain positivity of the cell average if the inflow boundary value and the source term are both positive, therefore the positivity-preserving framework in Zhang and Shu (2010) can be used to obtain a high order conservative positivity-preserving DG scheme. Unfortunately, in two-dimensions this is no longer the case. We show that the unmodulated DG solver based either on $$P^k$$ or $$Q^k$$ spaces (piecewise kth degree polynomials or piecewise tensor-product kth degree polynomials) could generate negative cell averages. We augment the DG space with additional functions so that the positivity of cell averages from the unmodulated DG solver can be restored, thereby leading to high order conservative positivity-preserving DG scheme based on these augmented DG spaces following the framework in Zhang and Shu (2010). Computational results are provided to demonstrate the good performance of our DG schemes.

We establish a rate of convergence for a semidiscrete operator splitting method applied to Hamilton--Jacobi equations with source terms. The method is based on sequentially solving a Hamilton--Jacobi equation and an ordinary differential equation. The Hamilton--Jacobi equation is solved exactly while the ordinary differential equation is solved exactly or by an explicit Euler method. We prove that the $L^{\infty}$ error associated with the operator splitting method is bounded by $\mathcal{O}(\Delta t)$, where $\Delta t$ is the splitting (or time) step. This error bound is an improvement over the existing $\mathcal{O}(\sqrt{\Delta t})$ bound due to Souganidis [Nonlinear Anal., 9 (1985), pp. 217--257]. In the one-dimensional case, we present a fully discrete splitting method based on an unconditionally stable front tracking method for homogeneous Hamilton--Jacobi equations. It is proved that this fully discrete splitting method possesses a linear convergence rate. Moreover, numerical results are presented to illustrate the theoretical convergence results.

Filtered schemes for Hamilton–Jacobi equations: A simple construction of convergent accurate difference schemes

With the increased integration density of modern radio frequency (RF) circuits and components, yield analysis and sensitivity study on design parameters have become ubiquitous in microwave engineering. To improve the electromagnetic simulation efficiency of such demands, a fast parameter sweep method based on the method of moment (MoM) in the planar layered medium is considered in this letter. By converting the sweeping of global parameters such as dielectric constants or substrate thickness to a matrix equation updating process in MoM, the overall complexity can be reduced by model order reduction (MOR) techniques. In addition, a representative circular patch antenna with a sectoral slot is studied and examined to demonstrate the effectiveness and accuracy of the proposed method.

High order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms

Computer implementation of an immersed boundary (IB) module inside a flow solver can be accomplished non-intrusively. However a versatile preprocessor is needed at the first place to extract the geometric information pertinent to the immersion of an arbitrarily complex geometry inside a Cartesian mesh. Geometric errors can negatively impact the correct implementation of the IB method as part of the solution algorithm. Additionally, the distance field from the geometry is needed to implement various turbulence models or flow initialization. Geometric processing stage for complex geometry have received less attention despite the popularity of the IB method. Our experience has shown that some of the procedures described in the literature have difficulties processing highly complex geometry or can be inflexible to implement reconstruction schemes for turbulent flows. To address these issues, we present a geometric preprocessor with a distance field solver. We constructed our procedure from computational geometry algorithms such as point-in-triangle and point-in-edge. The distance field solver uses the initial distance field at the immersed boundaries and propagates it to the rest of the domain by solving the Eikonal equation with the fast sweeping method. We demonstrate the versatility of our preprocessor for challenging test geometries from the computer graphics field, complex terrain and urban environments.

Traveltime computation is an effective way to simulate seismic wave propagation in isotropic and anisotropic media. It often requires the phase and group velocities along a ray direction. The phase and group velocities are not functions of the ray direction, but functions of the slowness direction, which motivates the need for efficient numerical approaches for computing the slowness direction. A generalized method was proposed to achieve this goal in 3D tilted transverse isotropic media, which involves solving a nonlinear equation of two unknowns, the inclination and azimuthal angles of the slowness direction, and is computationally demanding. To overcome the difficulty, the equation was recast by projecting it onto the local coordinates that align with the axis of symmetry of the medium. Under the local coordinate system, the nonlinear equation reduces to an equation of one unknown, the inclination angle, with the azimuthal angle obtained from the given ray direction. Hence, the nonlinear equation can be solved efficiently and its solutions can be used to calculate the phase and group velocities. As an application, the proposed method was used to compute traveltimes in 3D vertical transverse isotropic (VTI) media. The calculated group velocities are incorporated into the fast sweeping method to solve the original anisotropic Eikonal equation directly on uniform meshes. The feasibility of the proposed group velocity calculation method is evaluated in three 3D homogeneous anisotropic models, and the application on traveltime computation is tested on a 3D homogeneous VTI model and the British Petroleum anisotropic VTI model.

Computing the distance function to some surface or line is a problem that occurs very frequently. There are several ways of computing a relevant approximation of this function, using for example technique originating from the approximation of Hamilton Jacobi problems, or the fast sweeping method. Here we make a link with some elliptic problem and propose a very fast way to approximate the distance function.

Capturing stochastic variabilities in pedestrian flows: A dynamic continuum modeling approach

High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling