Solving the eikonal equation by the FSM method in Julia language

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.

Fast and accurate traveltime computation for quasi-P waves in anisotropic media is an essential ingredient of many seismic processing and interpretation applications such as Kirchhoff modeling and migration, microseismic source localization, and traveltime tomography. Fast-sweeping methods are widely used for solving the anisotropic eikonal equation due to their flexibility in solving general equations compared to the fast-marching method. However, it has been observed that fast sweeping can be much less efficient than fast marching for models with curved characteristics and practical grid sizes. By representing a tilted transversely isotropic (TTI) equation as a sequence of elliptically isotropic (EI) eikonal equations, we determine that the fast-marching algorithm can be used to compute fast and accurate traveltimes for TTI media. The tilt angle is absorbed into the description of the effective EI model; therefore, the adopted approach does not compromise on the solution accuracy. Through tests on benchmark synthetic models, we test our fast-marching algorithm and discover considerable improvement in accuracy by using factorization and a second-order finite-difference stencil. The adopted methodology opens the door to the possibility of using the fast-marching algorithm for a wider class of anisotropic eikonal equations.

A hybrid method for calculating seismic wave first-arrival traveltimes in two-dimensional models with an irregular surface

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.

A method towards robust and efficient medial axis transform (MAT) of arbitrary domains using distance solutions is presented. The distance field, d, is calculated by solving the hyperbolic-natured Eikonal (or Level Set) equation. The solution is obtained on Cartesian grids. Both the fast-marching method and fast-sweeping method are used to calculate d. Medial axis point clouds are then extracted based on the distance solution via a simple criteria: the Laplacian or the Hessian determinant of d. These point clouds in 2D-pixel and 3D-voxel space are further thinned to curves and surfaces through binary image thinning algorithms. This results in an overall hybrid approach. As an alternative to other methods, the current d −MAT procedure bypasses difficulties that are usually encountered by pure geometric methods (e.g. the Voronoi approach), especially in 3D, and provides better accuracy than pure thinning methods. It is also shown that the d −MAT approach provides the potential to sculpt/control the MAT form for specialized solution purposes. Various examples are given to demonstrate the current approach.

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.

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

We present StochasticPrograms.jl, a user-friendly and powerful open-source framework for stochastic programming written in the Julia language. The framework includes both modeling tools and structure-exploiting optimization algorithms. Stochastic programming models can be efficiently formulated using an expressive syntax, and models can be instantiated, inspected, and analyzed interactively. The framework scales seamlessly to distributed environments. Small instances of a model can be run locally to ensure correctness, whereas larger instances are automatically distributed in a memory-efficient way onto supercomputers or clouds and solved using parallel optimization algorithms. These structure-exploiting solvers are based on variations of the classical L-shaped, progressive-hedging, and quasi-gradient algorithms. We provide a concise mathematical background for the various tools and constructs available in the framework along with code listings exemplifying their usage. Both software innovations related to the implementation of the framework and algorithmic innovations related to the structured solvers are highlighted. We conclude by demonstrating strong scaling properties of the distributed algorithms on numerical benchmarks in a multinode setup. Summary of Contribution: This paper presents StochasticPrograms.jl, an open-source framework for stochastic programming implemented in the Julia programming language. The framework includes an expressive syntax for formulating stochastic programming models as well as versatile analysis tools and parallel optimization algorithms. The framework will prove useful to researchers, educators, and industrial users alike. Researchers will benefit from the readily extensible open-source framework, in which they can formulate complex stochastic models or quickly typeset and test novel optimization algorithms. Educators of stochastic programming will benefit from the clean and expressive syntax. Moreover, the framework supports analysis tools and stochastic programming constructs from classical theory and leading textbooks. We strongly believe that the StochasticPrograms.jl framework can reduce the barrier to entry for incoming practitioners of stochastic programming. Industrial practitioners can make use of StochasticPrograms.jl to rapidly formulate complex models, analyze small instances locally, and then run large-scale instances in production. In doing so, they get distributed capabilities for free without changing the code and access to well-tested state-of-the-art implementations of parallel structure-exploiting solvers. As the framework is open-source, anyone from these target audiences can contribute with new functionality to the framework. In conclusion, by providing both an intuitive interface for new users and an extensive development environment for expert users, StochasticPrograms.jl has strong potential to further the field of stochastic programming.

Summary High-frequency asymptotic methods, based on solving the eikonal equation, are widely used in many seismic applications including Kirchhoff migration and traveltime tomography. Finite-difference methods to solve the eikonal equation are computationally more efficient and attractive than ray tracing. But, finite-difference solution of the eikonal equation for a point-source contains inaccuracies due to source-singularity. Compared to the several proposed approaches to tackle source-singularity, factorization of the unknown traveltime is computationally efficient and simpler to implement. Recently, a factorization algorithm has been proposed to obtain clean first-order accuracy for tilted transversely isotropic (TTI) media. However, high-order accuracy of traveltimes is needed for computation of quantities that require traveltime derivatives, such as take-off angle and amplitudes. I propose an iterative fast sweeping algorithm to obtain high-order accuracy using factorization and a high-order finite-difference stencil. Numerical test shows improvements in accuracy of the TTI eikonal solution. This shows that once the source-singularity problem is tackled, high-order accurate solutions can be constructed easily. The method can be easily extended to media with lower anisotropic symmetries.

This paper develops a general method for solving reservoir simulation equations, which is based upon an implicit finite difference formulation of the partial differential equations of compositional thermal simulation, permitting the simulation of nearly all enhanced oil recovery techniques. The general method is used to derive 10 solution schemes, which include the well-known IMPES, Simultaneous Solution (SS), and Sequential Solution (SEQ) methods as sub-cases. The novelty of the general method lies in showing the interrelations of these methods, and also affording the possibility of deriving still other methods. The theoretical development of the general method is given in detail, and the special methods are derived in reference to the reservoir problems for which they were employed. The starting point is a system of nonlinear finite difference equations, typical of compositional thermal simulation. Newtonian iteration then leads to a system of equations. The partitioning of the left-hand side (containing the unknown vector) in a particular manner permits the derivation of 10 solution approaches, which are identified in detail. The proposed general method bridges the gap in the development of the reported solution methods, and also points to several new approaches. The procedure described can be utilized for a variety of simulation problems, with fully implicit source-sink terms (including heat loss, thermal cracking, etc.), for both the standard and multipoint difference schemes in multi dimensions.

A flux-limiting scheme for solution of the eikonal evolution equation

This chapter provides an overview of the mathematics of rays. It begins with a discussion of the theory of geometrical optics and how it can be formulated by means of the ray equations or the Hamilton-Jacobi equation. The two equations are of seemingly different types, but they are in fact equivalent. The ray equations are the characteristic equations of the Hamilton-Jacobi equation. This remark leads to the geometrical interpretation: the family of rays of geometrical optics is perpendicular to the wavefronts S = constant, if S denotes the appropriate solution of the Hamilton-Jacobi equation. The chapter considers the Hamilton-Jacobi theory in more detail, along with Hamilton's principle, ray differential geometry and the eikonal equation, and dispersion relations. It also presents the general solution of the linear wave equation before concluding with an analysis of the behavior of rays and waves in a slowly varying environment.

Reservoir simulators are highly complex computer programs and are often treated as “black boxes” that act on standardized input formats. In part, this is due to the closed-source nature of commercial offerings commonly used in industry, but even when the source-code is available, modification of a fully featured simulator can be a daunting task: A reservoir simulation problem, when fully specified, can have a great number of parameters and functional relationships that are tightly coupled together in most implementations. This is an obstacle to workflows that go beyond forecasting and into optimization, parameter fitting and sensitivity analysis where the solution quantity of interest should be paired with gradients taken with respect to forces and parameters. The monolithic design of typical reservoir simulation codes cannot be separated from the technical and mathematical choices made when discretizing the governing equations: The implicitness required to handle pressure changes and strong coupling between different solution variables means that the Jacobian of the governing equations must be implemented with respect to the primary variables: An error-prone process when done manually, and highly time-consuming when considering coupled systems. Automatic differentiation (AD) is a common cure to this problem that can generate Jacobians, but often comes at a severe cost in runtime performance. A more practical issue is that high computational performance for simulation problems relevant to modern geoenergy operations necessitates implementation in a compiled language while optimization and history matching workflows that integrate many data sources are most naturally expressed in high-level scripting languages. We describe the design of JutulDarcy, an open-source porous media simulator written in the Julia language. JutulDarcy is designed from the ground up with extensibility and computation of sensitivities in mind, starting from the assumption that high-performance automatic differentiation enables us to radically rethink the way simulators are designed. Highlights of this design include arbitrary execution graphs for the equation terms, in which dependence relationships can be re-configured on the fly, easy coupling of models, computation of sensitivities with respect to all model parameters through an adjoint method with a parameter declaration system, and an assembly process with performance that can exceed that of compiled simulators with hand-coded Jacobians. We validate the simulator on a number of standard test cases (immiscible, black-oil, CO2 and compositional flow) and demonstrate MPI parallel performance on multi-million cell models. In addition, we offer a frank assessment of the relative strengths and weaknesses of the Julia programming language for scientific computing.