Search

Conventional full-waveform inversion (FWI) using the least-squares norm ($L^2$) as a misfit function is known to suffer from cycle skipping. This increases the risk of computing a local rather than the global minimum of the misfit. In our previous work, we proposed the quadratic Wasserstein metric ($W_2$) as a new misfit function for FWI. The $W_2$ metric has been proved to have many ideal properties with regards to convexity and insensitivity to noise. When the observed and predicted seismic data are regarded as two density functions, the quadratic Wasserstein metric corresponds to the optimal cost of rearranging one density into the other, where the transportation cost is quadratic in distance. The difficulty of transforming seismic signals into nonnegative density functions is discussed. Unlike the $L^2$ norm, $W_2$ measures not only amplitude differences, but also global phase shifts, which helps to avoid cycle skipping issues. In this work, we build on our earlier method to cover more realistic high-resolution applications by embedding the $W_2$ technique into the framework of the adjoint-state method and applying it to seismic relevant 2D examples: the Camembert, the Marmousi, and the 2004 BP models. We propose a new way of using the $W_2$ metric trace-by-trace in FWI and compare it to global $W_2$ via the solution of the Monge-Amp\`ere equation. With corresponding adjoint source, the velocity model can be updated using the l-BFGS method. Numerical results show the effectiveness of $W_2$ for alleviating cycle skipping issues and sensitivity to noise. Both mathematical theory and numerical examples demonstrate that the quadratic Wasserstein metric is a good candidate for a misfit function in seismic inversion., Comment: 30 pages, 51 figures, 2nd version

We propose a new heterogeneous multiscale method (HMM) for computing the effective behavior of a class of highly oscillatory ordinary differential equations (ODEs). Without the need for identifying hidden slow variables, the proposed method is constructed based on the following ideas: a nonstandard splitting of the vector field (the right hand side of the ODEs); comparison of the solutions of the split equations; construction of effective paths in the state space whose projection onto the slow subspace has the correct dynamics; and a novel on-the-fly filtering technique for achieving a high order accuracy. Numerical examples are given.

The finite-difference time-domain method is one of the most popular for wave propagation in the time domain. One of its advantages is the use of a structured staggered grid, which makes it simple and efficient on modern computer architectures. A drawback, however, is the difficulty in approximating oblique boundaries, having to resort to staircase approximations. In many scattering problems this means that the grid resolution required to obtain an accurate solution is much higher than what is dictated by propagation in a homogeneous material. In this paper zero boundary data are considered, first for the velocity and then the pressure. These two forms of boundary conditions model perfectly rigid and pressure-release boundaries, respectively. A simple and efficient method to consistently model curved rigid boundaries in two dimensions was developed in Tornberg and Engquist [J. Comput. Phys. 227, 6922-6943 (2008)]. Here this treatment is generalized to three dimensions. Based on the approach of this method, a technique to model pressure-release surfaces with second order accuracy and without additional restriction on the timestep is also introduced. The structure of the standard method is preserved, making it easy to use in existing solvers. The effectiveness is demonstrated in several numerical tests.

This paper is concerned with preconditioning the stiffness matrix resulting from finite element discretizations of Maxwell's equations in the high frequency regime. The moving PML sweeping preconditioner, first introduced for the Helmholtz equation on a Cartesian finite difference grid, is generalized to an unstructured mesh with finite elements. The method dramatically reduces the number of GMRES iterations necessary for convergence, resulting in an almost linear complexity solver. Numerical examples including electromagnetic cloaking simulations are presented to demonstrate the efficiency of the proposed method.

This paper introduces a new sweeping preconditioner for the iterative solution of the variable coefficient Helmholtz equation in two and three dimensions. The algorithms follow the general structure of constructing an approximate $LDL^t$ factorization by eliminating the unknowns layer by layer starting from an absorbing layer or boundary condition. The central idea of this paper is to approximate the Schur complement matrices of the factorization using moving perfectly matched layers (PMLs) introduced in the interior of the domain. Applying each Schur complement matrix is equivalent to solving a quasi-1D problem with a banded LU factorization in the 2D case and to solving a quasi-2D problem with a multifrontal method in the 3D case. The resulting preconditioner has linear application cost and the preconditioned iterative solver converges in a number of iterations that is essentially indefinite of the number of unknowns or the frequency. Numerical results are presented in both two and three dimensions to demonstrate the efficiency of this new preconditioner., Comment: 25 pages

A new numerical method for approximating highly oscillatory wave fields as a superposition of Gaussian beams is presented. The method estimates the number of beams and their parameters automatically. This is achieved by an expectation-maximization algorithm that fits real, positive Gaussians to the energy of the highly oscillatory wave fields and its Fourier transform. Beam parameters are further refined by an optimization procedure that minimizes the difference between the Gaussian beam superposition and the highly oscillatory wave field in the energy norm.

In this paper, we develop and analyze a new class of discontinuous Galerkin (DG) methods for the acoustic wave equation in mixed form. Traditional mixed finite element (FE) methods produce energy conserving schemes, but these schemes are implicit, making the time-stepping inefficient. Standard DG methods give explicit schemes, but these approaches are typically dissipative or suboptimally convergent, depending on the choice of numerical fluxes. Our new method can be seen as a compromise between these two kinds of techniques, in the way that it is both explicit and energy conserving, locally and globally. Moreover, it can be seen as a generalized version of the Raviart-Thomas FE method and the finite volume method. Stability and convergence of the new method are rigorously analyzed, and we have shown that the method is optimally convergent. Furthermore, in order to apply the new method for unbounded domains, we apply our new method with the first order absorbing boundary condition. The stability of the resulting numerical scheme is analyzed.

A multiscale method for computing the effective behavior of a class of stiff and highly oscillatory ordinary differential equations (ODEs) is presented. The oscillations may be in resonance with on ...

A new set of consistent boundary conditions for Yee scheme approximations of wave equations in two space dimensions are developed and analyzed. We show how the classical staircase boundary conditions for hard reflections or, in the electromagnetic case, conducting surfaces in certain cases give O(1) errors. The proposed conditions keep the structure of the Yee scheme and are thus well suited for high performance computing. The higher accuracy is achieved by modifying the coefficients in the difference stencils near the boundary. This generalizes our earlier results with Gustafsson and Wahlund in one space dimension. We study stability and convergence and we present numerical examples.

Homogenization is an important mathematical framework for developing effective models of differential equations with oscillations. We include in the presentation techniques for deriving effective equations, a brief discussion on analysis of related limit processes and numerical methods that are based on homogenization principles. We concentrate on first- and second-order partial differential equations and present results concerning both periodic and random media for linear as well as nonlinear problems. In the numerical sections, we comment on computations of multi-scale problems in general and then focus on projection-based numerical homogenization and the heterogeneous multi-scale method.

A new approximation of the logarithmic derivative of the Hankel function is derived and applied to high-frequency wave scattering. We re-derive the on surface radiation condition (OSRC) approximations that are well suited for a Dirichlet boundary in acoustics. These correspond to the Engquist–Majda absorbing boundary conditions. Inverse OSRC approximations are derived and they are used for Neumann boundary conditions. We obtain an implicit OSRC condition, where we need to solve a tridiagonal system. The OSRC approximations are well suited for moderate wave numbers. The approximation of the logarithmic derivative is also used for deriving a generalized physical optics approximation, both for Dirichlet and Neumann boundary conditions. We have obtained similar approximations in electromagnetics, for a perfect electric conductor. Numerical computations are done for different objects in 2D and 3D and for different wave numbers. The improvement over the standard physical optics is verified.

Numerical approximation of wave propagation can be done very efficiently on uniform grids. The Yee scheme is a good example. A serious problem with uniform grids is the approximation of boundary co ...

Discretization of singular functions is an important component in many problems to which level set methods have been applied. We present two methods for constructing consistent approximations to Dirac delta measures concentrated on piecewise smooth curves or surfaces. Both methods are designed to be convenient for level set simulations and are introduced to replace the commonly used but inconsistent regularization technique that is solely based on a regularization parameter proportional to the mesh size. The first algorithm is based on a tensor product of regularized one-dimensional delta functions. It is independent of the irregularity relative to the grid. In the second method, the regularization is constructed from a one-dimensional regularization that is extended to multi-dimensions with a variable support depending on the orientation of the singularity relative to the computational grid. Convergence analysis and numerical results are given.

The heterogeneous multiscale methods (HMM) is a general framework for the numerical approximation of multiscale problems. It is here developed for ordinary differential equations containing differe ...

We will consider both explicit and implicit fully discrete finite volume schemes for solving three-dimensional Maxwell's equations with discontinuous physical coefficients on general polyhedral domains. Stability and convergence for both schemes are analyzed. We prove that the schemes are second order accurate in time. Both schemes are proved to be first order accurate in space for the Voronoi--Delaunay grids and second order accurate for nonuniform rectangular grids. We also derive explicit expressions for the dependence on the physical parameters in all estimates.

Singular terms in differential equations pose severe challenges for numerical approximations on regular grids. Regularization of the singularities is a very useful technique for their representation on the grid. We analyze such techniques for the practically preferred case of narrow support of the regularizations, extending our earlier results for wider support. The analysis also generalizes existing theory for one dimensional problems to multi-dimensions. New high order multi-dimensional techniques for differential equations and numerical quadrature are introduced based on the analysis and numerical results are presented. We also show that the common use of distance functions in level-set methods to extend one dimensional regularization to higher dimensions may produce O(1) errors.

Numerical simulation of high frequency acoustic, elastic or electro-magnetic wave propagation is important in many applications. Recently the traditional techniques of ray tracing based on geometrical optics have been augmented by numerical procedures based on partial differential equations. Direct simulations of solutions to the eikonal equation have been used in seismology, and lately approximations of the Liouville or Vlasov equation formulations of geometrical optics have generated impressive results. There are basically two techniques that follow from this latter approach: one is wave front methods and the other moment methods. We shall develop these methods in some detail after a brief review of more traditional algorithms for simulating high frequency wave propagation.

The rate of convergence for numerical methods approximating differential equations are often drastically reduced from lack of regularity in the solution. Typical examples are problems with singular source terms or discontinuous material coefficients. We shall discuss the technique of local regularization for handling these problems. New numerical methods are presented and analyzed and numerical examples are given. Some serious deficiencies in existing regularization methods are also pointed out.

Geometrical optics is a standard technique used for the approximation of high-frequency wave propagation. Computational methods based on partial differential equations instead of the traditional ray tracing have recently been applied to geometrical optics. These new methods have a number of advantages but typically exhibit difficulties with linear superposition of waves. In this paper we introduce a new partial differential technique based on the segment projection method in phase space. The superposition problem is perfectly resolved and so is the problem of computing amplitudes in the neighborhood of caustics. The computational complexity is of the same order as that of ray tracing. The new algorithm is described and a number of computational examples are given, including a simulation of waveguides.

We consider the inverse problem of determining the highly oscillatory coefficient $a^\epsilon$ in partial differential equations of the form $-\nabla\cdot (a^\epsilon\nabla u^\epsilon)+bu^\epsilon = f$ from given measurements of the solutions. Here, $\epsilon$ indicates the smallest characteristic wavelength in the problem ($0 0$, and exploration seismology, $b < 0$.

A systematic technique for the derivation of subgrid scale models in the numerical solution of partial differential equations is described. The technique is based on Haar wavelet projections of the discrete operator followed by a sparse approximation. As numerical testing suggests, the resulting numerical method will accurately represent subgrid scale phenomena on a coarse grid. Applications to numerical homogenization and wave propagation in materials with subgrid inhomogeneities are presented.

Finite difference approximations generically have ${\cal O}(1)$ pointwise errors close to a shock. We show that this local error may effect the smooth part of the solution such that only first order is achieved even for formally higher-order methods. Analytic and numerical examples of this form of accuracy are given. We also show that a converging method will have the formal order of accuracy in domains where no characteristics have passed through a shock.

In this paper we would like to point out some similarities between two artificial boundary conditions. One is the far field or absorbing boundary conditions for computations over unbounded domain. The other is the boundary conditions used at the boundary between subdomains in domain decomposition. We show some convergence result for the generalized Schwarz alternating method (GSAM), in which a convex combination of Dirichlet data and Neumann data is exchanged at the artificial boundary. We can see clearly how the mixed boundary condition and the relative size of overlap will affect the convergence rate. These results can be extended to more general coercive elliptic partial differential equations using the equivalence of elliptic operators. Numerically first- and second-order approximations of the Dirichlet-to-Neumann operator are constructed using local operators, where information tangential to the boundary is included. Some other possible extensions and applications are pointed out. Finally numerical results are presented.

A numerical homogenization procedure for elliptic differential equations is presented. It is based on wavelet decompositions of discrete operators in fine and coarse scale components followed by the elimination of the fine scale contributions. If the operator is in divergence form, this is preserved by the homogenization procedure. For periodic problems, results similar to classical effective coefficient theory are proved. The procedure can be applied to problems that are not cell-periodic.

Standard multigrid methods are not so effective for equations with highly oscillatory coefficients. New coarse grid operators based on homogenized operators are introduced to restore the fast convergence rate of multigrid methods. Finite difference approximations are used for the discretization of the equations. Convergence analysis is based on the homogenization theory. Proofs are given for a two-level multigrid method with the homogenized coarse grid operator for two classes of two-dimensional elliptic equations with Dirichlet boundary conditions.

We prove that for certain class of compactly supported C˜ initial data, smooth solutions of the unsteady Prandtl's equation blow up in nite time

New coarse grid operators are developed for elliptic problems with highly oscillatory coefficients. The new coarse grid operators are constructed directly based on the homogenized differential operators or hierarchically computed from the finest grid. A detailed description of this construction is provided. Numerical calculations for a two-dimensional elliptic model problem show that the homogenized form of the equations is very useful in the design of coarse grid operators for the multigrid method. A more realistic problem of heat conduction in a composite structure is also considered.

In this work we propose a new set of partial differential equations (PDEs) which can be seen as a generalization of the classical eikonal and transport equations, to allow for solutions with multiple phases. The traditional geometrical optics pair of equations suffer from the fact that the class of physically relevant solutions is limited. In particular, it does not include solutions with multiple phases, corresponding to crossing waves. Our objective has been to generalize these equations to accommodate solutions containing more than one phase. The new equations are based on the same high frequency approximation of the scalar wave equation as the eikonal and the transport equations. However, they also incorporate a finite superposition principle. The maximum allowed number of intersecting waves in the solution can be chosen arbitrarily, but a higher number means that a larger system of PDEs must be solved. The PDEs form a hyperbolic system of conservation laws with source terms.Although the equations are only weakly hyperbolic, and thus not well-posed in the strong sense, several examples show the viability of solving the equations numerically. The technique we use to capture multivalued solutions is based on a closure assumption for a system of equations representing the moments.

New numerical methods are derived for calculation of high frequency asymptotic expansion of the scalar wave equation. The nonlinear partial differential equations defining the terms in the expansio ...

A class was devised of fast wavelet based algorithms for linear evolution equations whose coefficients are time independent. The method draws on the work of Beylkin, Coifman, and Rokhlin which they applied to general Calderon-Zygmund type integral operators. A modification of their idea is applied to linear hyperbolic and parabolic equations, with spatially varying coefficients. A significant speedup over standard methods is obtained when applied to hyperbolic equations in one space dimension and parabolic equations in multidimensions.

Approximations of geometric optics type are commonly used in simulations of high frequency wave propagation. This form of technique fails when there is strong local variation in the wave speed on the scale of the wavelength or smaller. We propose a domain decomposition approach, coupling Gaussian beam methods where the wave speed is smooth with finite difference methods for the wave equations in domains with strong wave speed variation. In contrast to the standard domain decomposition algorithms, our finite difference domains follow the energy of the wave and change in time. A typical application in seismology presents a great simulation challenge involving the presence of irregularly located sharp inclusions on top of a smoothly varying background wave speed. These sharp inclusions are small compared to the domain size. Due to the scattering nature of the problem, these small inclusions will have a significant effect on the wave field. We present examples in two dimensions, but extensions to higher dimensions are straightforward.

In reflection seismology, seismic waves are recorded at the surface of the earth and migrated to image the subsurface. In mathematical terms, this is a boundary-value problem for the wave equation, where the boundary data are the recorded seismic wavefields. Gaussian beams, which are localized solutions of the wave equation, can be used for seismic migration. However, before Gaussian beams can be used for a solution of the boundary-value problem, the recorded data have to be represented in a compatible form. We present a new method for decomposing seismic data as a sum of Gaussian beams, based on transforming boundary data to a propagating wavefield. A simple synthetic seismic migration example illustrates our method and shows how different classes of Gaussian-beam parameters (such as curvature and extent) control the accuracy and sparseness of data representation.

TheHelmholtz equation describes wave propagation in the frequency domain and,as such, can be used for seismic imaging and fullwaveform inversion. We present two novel preconditioners for the efficientsolution of the Helmholtz equation in three dimensions. Both methodsfollow the general structure of constructing an approximate LDLt factorizationby eliminating the unknowns layer by layer starting from anabsorbing layer or boundary condition. In the first approach, werepresent the Schur complement matrices of the factorization in thehierarchical matrix framework. In the second approach, applying each Schurcomplement matrix is equivalent to solving a quasi-2D problem withthe multifrontal method. These preconditioners have linear application cost, andthe preconditioned iterative solvers converge in a number of iterationsthat is essentially independent of the number of unknowns orthe frequency. Numerical results on realistic 3D seismic models toconfirm the efficiency of these methods.

Difference approximations of hyperbolic partial differential equations with highly oscillatory coefficients and initial values are studied. Analysis of strong and weak convergence is carried out in the practically interesting case when the discretization step sizes are essentially independent of the oscillatory wave lengths. 01993 John Wiley & Sons, Inc.

The paper introduces the sweeping preconditioner, which is highly efficient for iterative solutions of the variable coefficient Helmholtz equation including very high frequency problems. The first central idea of this novel approach is to construct an approximate factorization of the discretized Helmholtz equation by sweeping the domain layer by layer, starting from an absorbing layer or boundary condition. Given this specific order of factorization, the second central idea of this approach is to represent the intermediate matrices in the hierarchical matrix framework. In two dimensions, both the construction and the application of the preconditioners are of linear complexity. The GMRES solver with the resulting preconditioner converges in an amazingly small number of iterations, which is essentially independent of the number of unknowns. This approach is also extended to the three dimensional case with some success. Numerical results are provided in both two and three dimensions to demonstrate the efficiency of this new approach., Comment: 34 pages

A triangle based total variation diminishing (TVD) scheme for the numerical approximation of hyperbolic conservation laws in two space dimensions is constructed. The novelty of the scheme lies in the nature of the preprocessing of the cell averaged data, which is accomplished via a nearest neighbor linear interpolation followed by a slope limiting procedures. Two such limiting procedures are suggested. The resulting method is considerably more simple than other triangle based non-oscillatory approximations which, like this scheme, approximate the flux up to second order accuracy. Numerical results for linear advection and Burgers' equation are presented.

This paper deals with numerical methods for high frequency wave scattering. It introduces a new hybrid technique that couples a directional fast multipole method for a subsection of a scattering surface to an asymptotic formulation over the rest of the scattering domain. The directional fast multipole method is new and highly efficient for the solution of the boundary integral formulation of a general scattering problem but it requires at least a few unknowns per wavelength on the boundary. The asymptotic method that was introduced by Bruno and collaborators requires much fewer unknowns. On the other hand the scattered field must have a simple structure. Hybridization of these two methods retains their best properties for the solution of the full problem. Numerical examples are given for the solution of the Helmholtz equation in two space dimensions.

discretized in time and space with step sizes At, £xx and Ay, t~ = n a t , xj = j£xx, Yk = kay. The numerical solution in the point (t~, xj, Yk) is denoted by uj~ and the entire solution at t ime t~ is u ~ and a function of j, k. We consider the solution of the above problem with finite difference methods. The problem is difficult because of the occurrence of discontinuities in the solution. One of the more widely used ways to derive a second order method with good shock capturing properties is to advance the solution with a first order method containing sufficient numerical damping ~, = a ( ~ ~)

A new class of computational far-field boundary conditions for hyperbolic partial differential equations was recently introduced by the authors. These boundary conditions combine properties of absorbing conditions for transient solutions and properties of far-field conditions for steady states. This paper analyses the properties of the wave equation coupled with these new boundary conditions: well-posedness, dissipativity and convergence in time.

We consider the wave equation in a medium with a rapidly varying speed of propagation. We construct a multiscale schemebased on the heterogeneous multiscale method, which can compute the correct coarse behavior of wave pulses traveling in themedium, at a computational cost essentially independent of the size of the small scale variations. This is verified by theoreticalresults and numerical examples. QC 20111117

Multiscale Discontinuous Galerkin Methods for Elliptic Problems with Multiple Scales.- Discrete Network Approximation for Highly-Packed Composites with Irregular Geometry in Three Dimensions.- Adaptive Monte Carlo Algorithms for Stopped Diffusion.- The Heterogeneous Multi-Scale Method for Homogenization Problems.- A Coarsening Multigrid Method for Flow in Heterogeneous Porous Media.- On the Modeling of Small Geometric Features in Computational Electromagnetics.- Coupling PDEs and SDEs: The Illustrative Example of the Multiscale Simulation of Viscoelastic Flows.- Adaptive Submodeling for Linear Elasticity Problems with Multiscale Geometric Features.- Adaptive Variational Multiscale Methods Based on A Posteriori Error Estimation: Duality Techniques for Elliptic Problems.- Multipole Solution of Electromagnetic Scattering Problems with Many, Parameter Dependent Incident Waves.- to Normal Multiresolution Approximation.- Combining the Gap-Tooth Scheme with Projective Integration: Patch Dynamics.- Multiple Time Scale Numerical Methods for the Inverted Pendulum Problem.- Multiscale Homogenization of the Navier-Stokes Equation.- Numerical Simulations of the Dynamics of Fiber Suspensions.

We shall consider systems of hyperbolic partial differential equations, $${{u}_{t}} + f{{(u)}_{x}} = 0.$$

Boundary conditions are derived for numerical wave simulation that minimize artificial reflections from the edges of the domain of computation. In this way acoustic and elastic wave propagation in a limited area can be efficiently used to describe physical behavior in an unbounded domain. The boundary conditions are based on paraxial approximations of the scalar and elastic wave equations. They are computationally inexpensive and simple to apply, and they reduce reflections over a wide range of incident angles.

Numerical approximations are studied for a large hyperbolic system coupled to a parabolic equation and a system of algebraic equations. The equations, which all are nonlinear, describe nonviscous compressible one-dimensional gas flow in a catalytic converter. Chemical reactions within the gas are included in the model. Well-posedness of the partial differential equations is analyzed together with stability of the numerical models. In particular an investigation is made of the effect of numerical dissipation and different boundary conditions. Numerical results are presented.

In recent years finite element methods have started to be applied to hyperbolic equations. Since modern finite element and finite difference methods for hyperbolic equations look very much alike, new results in the analysis of difference methods are also applicable to element methods. We shall discuss propagation of sharp signals, problems with different time scales and the effect of boundaries on stability and accuracy.

In practical calculations, it is often essential to introduce artificial boundaries to limit the area of computation. Here we develop a systematic method for obtaining a hierarchy of local boundary conditions at these artificial boundaries. These boundary conditions not only guarantee stable difference approximations but also minimize the (unphysical) artificial reflections which occur at the boundaries.