Search

We propose a novel Babich-like ansatz consisting of an infinite series of dyadic coefficients (three-by-three matrices) and spherical Hankel functions for solving point-source Maxwell's equations in an inhomogeneous medium so as to produce the so-called dyadic Green's function. Using properties of spherical Hankel functions, we derive governing equations for the unknown asymptotics of the ansatz including the traveltime function and dyadic coeffcients. By proposing matching conditions at the point source, we rigorously derive asymptotic behaviors of these geometrical-optics ingredients near the source so that their initial data at the source point are well-defined. To verify the feasibility of the proposed ansatz, we truncate the ansatz to keep only the first two terms, and we further develop partial differential equation-based Eulerian approaches to compute the resulting asymptotic solutions. Since the system of governing equations for each dyadic coeffcient is strongly coupled, we introduce auxiliary variables to transform these strongly coupled systems into decoupled scalar equations. Furthermore, we develop high-order Lax-Friedrichs weighted essentially nonoscillatory schemes for computing these auxiliary variables so that the Green's function can be constructed. Numerical examples demonstrate that our new ansatz yields a uniform asymptotic solution in the region of space containing a point source but no other caustics. [ABSTRACT FROM AUTHOR]

In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation so that Maxwell's equations may be viewed as an evolution equation in one of the spatial directions. With such applications in mind, we propose a new Eulerian geometrical-optics method, dubbed the fast Huygens sweeping method, for computing Green's functions of Maxwell's equations in inhomogeneous media in the high-frequency regime and in the presence of caustics. The first novelty of the fast Huygens sweeping method is that a new dyadic-tensor-type geometrical-optics ansatz is proposed for Green's functions which is able to utilize some unique features of Maxwell's equations. The second novelty is that the Huygens{Kirchhoff secondary source principle is used to integrate many locally valid asymptotic solutions to yield a globally valid asymptotic solution so that caustics associated with the usual geometrical-optics ansatz can be treated automatically. The third novelty is that a buttery algorithm is adapted to carry out the matrix-vector products induced by the Huygens{Kirchhoff integration in O(N logN) operations, where N is the total number of mesh points, and the proportionality constant depends on the desired accuracy and is independent of the frequency parameter. To reduce the storage of the resulting traveltime and amplitude tables, we compress each table into a linear combination of tensor-product based multivariate Chebyshev polynomials so that the information of each table is encoded into a small number of Chebyshev coeffcients. The new method enjoys the following desired features: (1) it precomputes a set of local traveltime and amplitude tables; (2) it automatically takes care of caustics; (3) it constructs Green's functions of Maxwell's equations for arbitrary frequencies and for many point sources; (4) for a specified number of points per wavelength it constructs each Green's function in nearly optimal complexity O(N logN) in terms of the total number of mesh points N, where the prefactor of the complexity depends only on the specified accuracy and is independent of the frequency parameter. Three-dimensional numerical experiments are presented to demonstrate the performance and accuracy of the new method. [ABSTRACT FROM AUTHOR]

We have developed a level-set method for the inverse gravimetry problem of imaging salt structures with density contrast reversal. Under such a circumstance, a part of the salt structure contributes two completely opposite anomalies that counteract with each other, making it unobservable to the gravity data. As a consequence, this amplifies the inherent nonuniqueness of the inverse gravimetry problem so that it is much more challenging to recover the whole salt structure from the gravity data. To alleviate the severe nonuniqueness, it is reasonable to assume that the density contrast between the salt structure and the surrounding sedimentary host depends upon the depth only and is known a priori. Consequently, the original inverse gravity problem reduces to a domain inverse problem, where the supporting domain of the salt body becomes the only unknown. We have used a level-set function to parametrize the boundary of the salt body so that we reformulated the domain inverse problem into a nonlinear optimization problem for the level-set function, which was further solved for by a gradient descent method. Both 2D and 3D experiments on the SEG/EAGE salt model were carried out to demonstrate the effectiveness and efficiency of the new method. The algorithm was able to recover dipping flanks of the salt model, and it only took 40 min in a 2.5 GHz CPU to invert for a 3D model of 97,000 unknowns. [ABSTRACT FROM AUTHOR]

Wave propagation in an isotropic acoustic medium occupied by a moving fluid is governed by an anisotropic eikonal equation. Since this anisotropic eikonal equation is associated with an inhomogeneous Hamiltonian, most of existing anisotropic eikonal solvers are either inapplicable or of unpredictable behavior in convergence. Realizing that this anisotropic eikonal equation is defined by a sum of two well-understood first-order differential operators, we propose novel operator-splitting based fast sweeping methods to solve this generalized eikonal equation. We develop various operator-splitting methods relying on the Peaceman-Rachford scheme, the Douglas- Rachford scheme, the θ-scheme, and the regularized θ-scheme. After applying the operator-splitting strategy, each splitting step corresponds to a much simpler Hamilton-Jacobi equation so that we can apply the Lax-Friedrichs sweeping method to solve these splitted equations efficiently and easily. Two- and three-dimensional examples demonstrate the performances and efficiency of the proposed approaches. [ABSTRACT FROM AUTHOR]

Accurate simulation of seismic waves is of critical importance in a variety of geophysical applications. Based on recent works on staggered discontinuous Galerkin methods, we have developed a new method for the simulations of seismic waves, which has energy conservation and extremely low grid dispersion, so that it naturally provided accurate numerical simulations of wave propagation useful for geophysical applications and was a generalization of classical staggered-grid finite-difference methods. Moreover, it could handle with ease irregular surface topography and discontinuities in the subsurface models. Our new method discretized the velocity and the stress tensor on this staggered grid, with continuity imposed on different parts of the mesh. The symmetry of the stress tensor was enforced by the Lagrange multiplier technique. The resulting method was an explicit scheme, requiring the solutions of a block diagonal system and a local saddle point system in each time step, and it was, therefore, very efficient. To tailor our scheme to Rayleigh waves, we developed a mortar formulation of our method. Specifically, a fine mesh was used near the free surface and a coarse mesh was used in the rest of the domain. The two meshes were in general not matching, and the continuity of the velocity at the interface was enforced by a Lagrange multiplier. The resulting method was also efficient in time marching. We also developed a stability analysis of the scheme and an explicit bound for the time step size. In addition, we evaluated some numerical results and found that our method was able to preserve the wave energy and accurately computed the Rayleigh waves. Moreover, the mortar formulation gave a significant speed up compared with the use of a uniform fine mesh, and provided an efficient tool for the simulation of Rayleigh waves. [ABSTRACT FROM AUTHOR]

We propose a new methodology for carrying out eikonal based traveltime tomography arising from important applications such as seismic imaging and medical imaging. The new method formulates the traveltime tomography problem as a variational problem for a certain cost functional explicitly with respect to both traveltime and sound speed. Furthermore, the cost functional is penalized to enforce the nonlinear equality constraint associated with the underlying eikonal equation, biharmonically regularized with respect to traveltime, and harmonically regularized with respect to sound speed. To overcome the difficulty associated with the inherent nonlinearity of the eikonal equation, the Euler--Lagrange equation of the penalized-regularized variational problem is reformulated into an equivalent, mixed optimality system. This mixed system is associated with an initial value problem which is solved by an operator-splitting based solution method, and the splitting approach effectively reduces the optimality system into three nonlinear subproblems and three linear subproblems. Moreover, the nonlinear subproblems can be solved pointwise, while the linear subproblems can be reduced to linear second-order elliptic problems. Numerical experiments show that the new method can carry out traveltime tomography successfully and recover sound speeds efficiently. [ABSTRACT FROM AUTHOR]

We have developed a local level-set method for inverting 3D gravity-gradient data. To alleviate the inherent non-uniqueness of the inverse gradiometry problem, we assumed that a homogeneous density contrast distribution with the value of the density contrast specified a priori was supported on an unknown bounded domain D so that we may convert the original inverse problem into a domain inverse problem. Because the unknown domain D may take a variety of shapes, we parametrized the domain D by a level-set function implicitly so that the domain inverse problem was reduced to a nonlinear optimization problem for the level-set function. Because the convergence of the level-set algorithm relied heavily on initializing the level-set function to enclose the gravity center of a source body, we applied a weighted L¹-regularization method to locate such a gravity center so that the level-set function can be properly initialized. To rapidly compute the gradient of the nonlinear functional arising in the level-set formulation, we made use of the fact that the Laplacian kernel in the gravity force relation decayed rapidly off the diagonal so that matrix-vector multiplications for evaluating the gradient can be accelerated significantly. We conducted extensive numerical experiments to test the performance and effectiveness of the new method. [ABSTRACT FROM AUTHOR]

Motivated by recent theoretical results obtained by the third author for the identification problem arising in single-photon emission computerized tomography (SPECT), we propose an adjoint state method for recovering both the source and the attenuation in the attenuated X-ray transform. Our starting point is the transport-equation characterization of the attenuated X-ray transform, and we apply efficient fast sweeping methods to solve static transport equations and adjoint state equations. Numerous examples are presented to demonstrate various features of the identification problem, such as uniqueness and nonuniqueness, stability and instability, and recovery of the wave front set. [ABSTRACT FROM AUTHOR]

The solution for the eikonal equation with a point-source condition has an upwind singularity at the source point as the eikonal solution behaves like a distance function at and near the source. As such, the eikonal function is not differentiable at the source so that all formally high-order numerical schemes for the eikonal equation yield first-order convergence and relatively large errors. Therefore, it is a longstanding challenge in computational geometrical optics how to compute a uniformly high-order accurate solution for the point-source eikonal equation in a global domain. In this paper, assuming that both the squared slowness and the squared eikonal are analytic near the source, we propose high-order factorization based high-order hybrid fast sweeping methods for point-source eikonal equations to compute just such solutions. Observing that the squared eikonal is differentiable at the source, we propose to factorize the eikonal into two multiplicative or additive factors, one of which is specified to approximate the eikonal up to arbitrary order of accuracy near the source, and the other of which serves as a higher-order correction term. This decomposition is achieved by using the eikonal equation and applying power series expansions to both the squared eikonal and the squared slowness function. We develop recursive formulas to compute the approximate eikonal up to arbitrary order of accuracy near the source. Furthermore, these approximations enable us to decompose the eikonal into two factors, either multiplicatively or additively, so that we can design two new types of hybrid, high-order fast sweeping schemes for the point-source eikonal equation. We also show that the first-order hybrid fast sweeping methods are monotone and consistent so that they are convergent in computing viscosity solutions. Two- and three-dimensional numerical examples demonstrate that a hybrid pth order fast sweeping method yields desired, uniform, clean pth order convergence in a global domain by using a pth order factorization. [ABSTRACT FROM AUTHOR]

In this paper, we obtain the first local a posteriori error estimate for time-dependent Hamilton-Jacobi equations. Given an arbitrary domain Ω and a time T, the estimate gives an upper bound for the L∞-norm in Ω at time T of the difference between the viscosity solution u and any continuous function v in terms of the initial error in the domain of dependence and in terms of the (shifted) residual of v in the union of all the cones of dependence with vertices in Ω. The estimate holds for general Hamiltonians and any space dimension. It is thus an ideal tool for devising adaptive algorithms with rigorous error control for time-dependent Hamilton-Jacobi equations. This result is an extension to the global a posteriori error estimate obtained by S. Albert, B. Cockburn, D. French, and T. Peterson in A posteriori error estimates for general numerical methods for Hamilton-Jacobi equations. Part II: The time-dependent case, Finite Volumes for Complex Applications, vol. III, June 2002, pp. 17-24. Numerical experiments investigating the sharpness of the a posteriori error estimates are given. [ABSTRACT FROM AUTHOR]

In the geometrical-optics approximation for the Helmholtz equation with a point source, traveltimes and amplitudes have upwind singularities at the point source. Hence, both first-order and higher-order finite-difference solvers exhibit formally at most first-order convergence and relatively large errors. Such singularities can be factored out by factorizing traveltimes and amplitudes, where one factor is specified to capture the corresponding source singularity and the other factor is an unknown function smooth near the source. The resulting underlying unknown functions satisfy factored eikonal and transport equations, respectively. A third-order Lax-Friedrichs scheme is designed to compute the underlying functions. Thus, highly accurate first-arrival traveltimes and reliable amplitudes can be computed. Furthermore, asymptotic wavefields are constructed using computed traveltimes and amplitudes in the WKBJ form. Two-dimensional and 3D examples demonstrate the performance of the proposed algorithms, and the constructed WKBJ Green's functions are in good agreement with direct solutions of the Helmholtz equation before caustics occur. [ABSTRACT FROM AUTHOR]

We prove optimal convergence rates for the approximation provided by the original discontinuous Galerkin method for the transport-reaction problem. This is achieved in any dimension on meshes related in a suitable way to the possibly variable velocity carrying out the transport. Thus, if the method uses polynomials of degree k, the L2-norm of the error is of order k+1. Moreover, we also show that, by means of an element-by-element postprocessing, a new approximate flux can be obtained which superconverges with order k + 1. [ABSTRACT FROM AUTHOR]

We design an Eulerian Gaussian beam summation method for solving Helmholtz equations in the high-frequency regime. The traditional Gaussian beam summation method is based on Lagrangian ray tracing and local ray-centered coordinates. We propose a new Eulerian formulation of Gaussian beam theory which adopts global Cartesian coordinates, level sets, and Liouville equations, yielding uniformly distributed Eulerian traveltimes and amplitudes in phase space simultaneously for multiple sources. The time harmonic wavefield can be constructed by summing up Gaussian beams with ingredients provided by the new Eulerian formulation. The conventional Gaussian beam summation method can be derived from the proposed method. There are three advantages of the new method: (1) We have uniform resolution of ray distribution. (2) We can obtain wavefields excited at different sources by varying only source locations in the summation formula. (3) We can obtain wavefields excited at different frequencies by varying only frequencies in the summation formula. Numerical experiments indicate that the Gaussian beam summation method yields accurate asymptotic wavefields even at caustics. The new method may be used for seismic modeling and migration. [ABSTRACT FROM AUTHOR]

A new numerical procedure coupling the level-set method with the moving-mesh method to simulate subcooled nucleate pool boiling is proposed. Numerical test problems have validated this new method. The simulation of bubble dynamics during nucleate boiling under liquid subcooling shows that this novel adaptive method is more accurate in determining interfacial heat transfer than a computational method based on uniform grids with the same number of mesh points. [ABSTRACT FROM AUTHOR]

The original fast sweeping method, which is an efficient iterative method for stationary Hamilton--Jacobi equations, relies on natural ordering provided by a rectangular mesh. We propose novel ordering strategies so that the fast sweeping method can be extended efficiently and easily to any unstructured mesh. To that end we introduce multiple reference points and order all the nodes according to their lp-metrics to those reference points. We show that these orderings satisfy the two most important properties underlying the fast sweeping method: (1) these orderings can cover all directions of information propagating efficiently; (2) any characteristic can be decomposed into a finite number of pieces and each piece can be covered by one of the orderings. We prove the convergence of the new algorithm. The computational complexity of the algorithm is nearly optimal in the sense that the total computational cost consists of O(M) flops for iteration steps and O(MlogM) flops for sorting at the predetermined initialization step which can be efficiently optimized by adopting a linear time sorting method, where M is the total number of mesh points. Extensive numerical examples demonstrate that the new algorithm converges in a finite number of iterations independent of mesh size. [ABSTRACT FROM AUTHOR]

We propose a local level set method for constructing the geometrical optics term in the paraxial formulation for the high frequency asymptotics of two-dimensional (2-D) acoustic wave equations. The geometrical optics term consists of two multivalued functions: a travel-time function satisfying the eikonal equation locally and an amplitude function solving a transport equation locally. The multivalued travel-times are obtained by solving a level set equation and a travel-time equation with a forcing term. The multivalued amplitudes are computed by a new Eulerian formula based on the gradients of travel-times and takeoff angles. As a byproduct the method is also able to capture the caustic locations. The proposed Eulerian method has complexity of O(N2LogN), rather than O(N4) as typically seen in the Lagrangian ray tracing method. Several examples including the well-known Marmousi synthetic model illustrate the accuracy and efficiency of the Eulerian method. [ABSTRACT FROM AUTHOR]

Traveltime, or geodesic distance, is locally the solution of the eikonal equation of geometric optics. However traveltime between sufficiently distant points is generically multivalued. Finite difference eikonal solvers approximate only the viscosity solution, which is the smallest value of the (multivalued) traveltime ("first arrival time"). The slowness matching method stitches together local single-valued eikonal solutions, approximated by a finite difference eikonal solver, to approximate all values of the traveltime. In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation, so that the eikonal equation may be viewed as an evolution equation in one of the spatial directions. This paraxial assumption simplifies both the efficient computation of local traveltime fields and their combination into global multivalued traveltime fields via the slowness matching algorithm. The cost of slowness matching is on the same order as that of a finite difference solver used to compute the viscosity solution, when traveltimes from many point sources are required as is typical in seismic applications. Adaptive gridding near the source point and a formally third order scheme for the paraxial eikonal combine to give second order convergence of the traveltime branches. [ABSTRACT FROM AUTHOR]

In this work, an adaptive strategy for the phase space method for traveltime tomography (Chung et al 2007 Inverse Problems 23 309-29) is developed. The method first uses those geodesics/rays that produce smaller mismatch with the measurements and continues on in the spirit of layer stripping without defining the layers explicitly. The adaptive approach improves stability, efficiency and accuracy. We then extend our method to reflection traveltime tomography by incorporating broken geodesics/rays for which a jump condition has to be imposed at the broken point for the geodesic flow. In particular, we show that our method can distinguish non-broken and broken geodesics in the measurement and utilize them accordingly in reflection traveltime tomography. We demonstrate that our method can recover the convex hull (with respect to the underlying metric) of unknown obstacles as well as the metric outside the convex hull. [ABSTRACT FROM AUTHOR]