Computation Of Wave Fields Research Articles

Closing the performance gap between an iterative frequency-domain solver and an explicit time-domain scheme for 3D migration on parallel architectures

Three-dimensional reverse-time migration with the constant-density acoustic wave equation requires an efficient numerical scheme for the computation of wavefields. An explicit finite-difference scheme in the time domain is a common choice. However, it requires a significant amount of disk space for the imaging condition. The frequency-domain approach simplifies the correlation of the source and receiver wavefields, but requires the solution of a large sparse linear system of equations. For the latter, we use an iterative Krylov solver based on a shifted Laplace multigrid preconditioner with matrix-dependent prolongation. The question is whether migration in the frequency domain can compete with a time-domain implementation when both are performed on a parallel architecture. Both methods are naturally parallel over shots, but the frequency-domain method is also parallel over frequencies. If we have a sufficiently large number of compute nodes, we can compute the result for each frequency in parallel and the required time is dominated by the number of iterations for the highest frequency. As a parallel architecture, we consider a commodity hardware cluster that consists of multicore central processing units (CPUs), each of them connected to two graphics processing units (GPUs). Here, GPUs are used as accelerators and not as an independent compute node. The parallel implementation of the 3D migration in frequency domain is compared to a time-domain implementation. We optimize the throughput of the latter with dynamic load balancing, asynchronous I/O, and compression of snapshots. Because the frequency-domain solver uses matrix-dependent prolongation, the coarse-grid operators require more storage than available on GPUs for problems of realistic size. Due to data transfer, there is no significant speedup using GPU-accelerators. Therefore, we consider an implementation on CPUs only. Nevertheless, with the parallelization over shots and frequencies, this approach could compete with the time-domain implementation on multiple GPUs.

Stability for Imposing Absorbing Boundary Conditions in the Finite Element Simulation of Acoustic Wave Propagation

It is well-known that artificial boundary conditions are crucial for the efficient and accurate computations of wavefields on unbounded domains. In this paper, we investigate stability analysis for the wave equation coupled with the first and the second order absorbing boundary conditions. The computational scheme is also developed. The approach allows the absorbing boundary conditions to be naturally imposed, which makes it easier for us to construct high order schemes for the absorbing boundary conditions. A thirdorder Lagrange finite element method with mass lumping is applied to obtain the spatial discretization of the wave equation. The resulting scheme is stable and is very efficient since no matrix inversion is needed at each time step. Moreover, we have shown both abstract and explicit conditional stability results for the fully-discrete schemes. The results are helpful for designing computational parameters in computations. Numerical computations are illustrated to show the efficiency and accuracy of our method. In particular, essentially no boundary reflection is seen at the artificial boundaries. Mathematics subject classification: 35L05, 35L20, 65M06, 65M12, 65M60.

Lithospheric expression of geological units in central and eastern North America from full waveform tomography

The EarthScope TA deployment has provided dense array coverage throughout the continental US and with it, the opportunity for high resolution 3D seismic velocity imaging of both lithosphere and asthenosphere in the continent. Building upon our previous long-period waveform tomographic modeling in North America, we present a higher resolution 3D isotropic and radially anisotropic shear wave velocity model of the North American lithospheric mantle, constructed tomographically using the spectral element method for wavefield computations and waveform data down to 40 s period. The new model exhibits pronounced spatial correlation between lateral variations in seismic velocity and anisotropy and major tectonic units as defined from surface geology. In the center of the continent, the North American craton exhibits uniformly thick lithosphere down to 200–250 km, while major tectonic sutures of Proterozoic age visible in the surface geology extend down to 100–150 km as relatively narrow zones of distinct radial anisotropy, with Vsv>Vsh. Notably, the upper mantle low velocity zone is present everywhere under the craton between 200 and 300 km depth. East of the continental rift margin, the lithosphere is broken up into a series of large, somewhat thinner (150 km) high velocity blocks, which extend laterally 200–300 km offshore into the Atlantic Ocean. Between the craton and these deep-rooted blocks, we find a prominent narrow band of low velocities that roughly follows the southern and eastern Laurentia rift margin and extends into New England. We suggest that the lithosphere along this band of low velocities may be thinned due to the combined effects of repeated rifting processes and northward extension of the hotspot related Bermuda low-velocity channel across the New England region. We propose that the deep rooted high velocity blocks east of the Laurentia margin represent the Proterozoic Gondwanian terranes of pan-African affinity, which were captured during the Rodinia formation but left behind after the opening of the Atlantic Ocean. Our results suggest that recurring episodes of tectonic events that are well exposed at the surface also leave persistent scars in the continental lithosphere mantle, marked by isotropic and radially anisotropic velocity anomalies that reach as deep as 100–150 km.

A multilevel fast multiple method for computing the propagation of multiply scattered 2.5-D teleseismic surface waves underneath a linear or quasi-linear seismic station array

We introduce an algorithm for the forward modeling of multiple scattering of teleseismic surface waves where the underground structure beneath the linear (or quasi-linear) seismic station array is assumed two-dimensional and the teleseismic surface waves may be approached from an arbitrary direction. The current algorithm is two-and-half-dimensional since even though the structure is two-dimensional, the displacements are three-dimensional because of scattering outside the sagittal plane. The total displacement field is computed by applying a convolutional type integral equation for which the Green’s function is obtained by employing the normal mode theory in one-dimensional reference structure. When the number of data points gets large, the wavefield computations become awkward. In order to alleviate the computational burden, we employ the proficient multilevel fast multipole method where the central processing unit (CPU) time increases logarithmically with the size of the model. Key words: Green’s function, heterogeneity, integral equation, multiple scattering, multipole expansions, surface waves

Modeling of Lorentz forces and radiated wave fields for bulk wave electromagnetic acoustic transducers

Currently, the finite element method (FEM) and analytical calculation are widely employed for the modeling of electromagnetic acoustic transducers (EMATs). However, it takes long time for finite element calculation. Previous analytical models for bulk wave EMATs are generally considered separately and incompletely, and expressions of radiated wave fields contain infinite integrations and multiple singular points, which result in complex numerical computation. A complete model containing the Lorentz force and radiated wave field calculation for the EMAT with a spiral coil and a NdFeB permanent magnet is established. By introducing a current loop instead of the permanent magnet and adopting the truncated region eigenfunction expansion (TREE) method, the distributions of static and dynamic magnetic fields and their generated Lorentz forces are calculated. A series expansion method is proposed for the computation of radiated wave fields, which replaces the integration by series operation and avoids the solutions of singular points effectively. The Lorentz forces and radiated wave fields of a typical transducer are computed. The validity of the model is verified by FEM and experiments. Their good agreements verify the accuracy and validity of the model.

An Educational Resource for Visualizing the Global Seismic Wave Field

The technology currently available in university classrooms often includes interactive whiteboards, high‐definition displays and projectors; furthermore, many students are equipped with tablet computers. Nevertheless, teaching material that can take advantage of these technologies is time consuming to produce and there is need for additional available material for courses in geophysics. The use of dynamic materials, such as animations and movies in the classroom has been shown to have the ability to enhance learning (e.g., see Mayer and Moreno, 2002, for a review). However, traditional methods of teaching about the global seismic wave field typically involve static imagery. Usually this has been achieved through a combination of drawing ray paths and wave fronts at distinct timesteps. Thus, students have had to build the mental connections between the underlying wave propagation and the static imagery being displayed. But the recent availability of cheap computer clusters has made numerical computation of the seismic wave field in realistic Earth models readily accessible. As a result, animating seismic‐wave propagation has become relatively commonplace. Some of the first global animations of seismic‐wave propagation were produced by summing torsional mode free oscillations in the mantle to show the SH wave field at periods down to 12 s (Wysession and Shore, 1994). This set of animations has been a valuable educational resource and has inspired us to create a new collection. Here, we present a new series of global seismic‐wave animations that includes both the SH and P / SV wave field. In addition, recent developments to the TauP Toolkit (Crotwell et al. , 1999) have made it possible to animate seismic ray paths and wavefronts. Thus, our collection combines the seismic wave field with dynamically drawn ray paths and/or wavefronts. In what follows, we describe a few of the animations that are available on the global scale and …

Mathematical simulation of acoustic wave refraction near a caustic

Issues related to the computation of wave fields in an acoustic medium near caustics are considered. A boundary condition on a caustic is established, and the Green's function of a boundary value problem for the general case of a varying speed of sound is constructed. For this purpose, an auxiliary Goursat problem is considered and a system of its particular solutions is constructed using hyper-geometric functions. A Volterra integral equation for the Green's function is obtained, and an algorithm for its expansion with respect to smoothness is described. A finite difference scheme approximating the solution of the differential problem with an unbounded coefficient is proposed. Numerical results are presented.

A fast parallel algorithm for solving block-tridiagonal systems of linear equations including the domain decomposition method

In this study, we develop a new parallel algorithm for solving systems of linear algebraic equations with the same block-tridiagonal matrix but with different right-hand sides. The method is a generalization of the parallel dichotomy algorithm for solving systems of linear equations with tridiagonal matrices \cite{terekhov:Dichotomy}. Using this approach, we propose a parallel realization of the domain decomposition method (\mbox{the Schur} complement method). The calculation of acoustic wave fields using the spectral-difference technique improves the efficiency of the parallel algorithms. A near-linear dependence of the speedup with the number of processors is attained using both several and several thousands of processors. This study is innovative because the parallel algorithm developed for solving block-tridiagonal systems of equations is an effective and simple set of procedures for solving engineering tasks on a supercomputer.

Reliability assessment of different plate theories for elastic wave propagation analysis in functionally graded plates

The importance of elastic wave propagation problem in plates arises from the application of ultrasonic elastic waves in non-destructive evaluation of plate-like structures. However, precise study and analysis of acoustic guided waves especially in non-homogeneous waveguides such as functionally graded plates are so complicated that exact elastodynamic methods are rarely employed in practical applications. Thus, the simple approximate plate theories have attracted much interest for the calculation of wave fields in FGM plates. Therefore, in the current research, the classical plate theory (CPT), first-order shear deformation theory (FSDT) and third-order shear deformation theory (TSDT) are used to obtain the transient responses of flexural waves in FGM plates subjected to transverse impulsive loadings. Moreover, comparing the results with those based on a well recognized hybrid numerical method (HNM), we examine the accuracy of the plate theories for several plates of various thicknesses under excitations of different frequencies. The material properties of the plate are assumed to vary across the plate thickness according to a simple power-law distribution in terms of volume fractions of constituents. In all analyses, spatial Fourier transform together with modal analysis are applied to compute displacement responses of the plates. A comparison of the results demonstrates the reliability ranges of the approximate plate theories for elastic wave propagation analysis in FGM plates. Furthermore, based on various examples, it is shown that whenever the plate theories are used within the appropriate ranges of plate thickness and frequency content, solution process in wave number-time domain based on modal analysis approach is not only sufficient but also efficient for finding the transient waveforms in FGM plates.

Wave Trains Associated with a Cascade of Bifurcations of Space-Time Caustics over Elongated Underwater Banks

We study the behavior of linear nonstationary shallow water waves generated by an instantaneous localized source as they propagate over and become trapped by elongated underwater banks or ridges. To find the solutions of the corresponding equations, we use an earlier-developed asymptotic approach based on a generalization of Maslov’s canonical operator, which provides a relatively simple and efficient analytic-numerical algorithm for the wave field computation. An analysis of simple examples (where the bank and source shapes are given by certain elementary functions) shows that the appearance and dynamics of trapped wave trains is closely related to a cascade of bifurcations of space-time caustics, the bifurcation parameter being the bank length-to-width ratio.

Site characterization using Gauss–Newton inversion of 2-D full seismic waveform in the time domain

Site characterization for design of deep foundations is very crucial, as unanticipated site conditions still represent significant problems and disputes occur during construction. Traditional surface-based geophysical methods, which use wave velocity dispersion or first-arrival times, have been widely used recently to assess spatial variation; however they cannot well characterize reverse profiles or buried low-velocity zones. For better characterization of these challenging site conditions, a full waveform inversion based on Gauss–Newton method is presented. The inversion scheme is based on a finite-difference solution of the 2-D elastic wave equation in the time domain. The strength of this approach is the ability to generate all possible wave types of seismic wavefields that are then compared with observed data to infer complex subsurface properties. Virtual sources and reciprocity of wavefields are used for calculation of partial derivative wavefields to reduce computer time. Cross convolution between observed and estimated wavefields are also employed to allow the technique to be independent of the source signatures. The capability of the presented technique is tested with both synthetic and real experimental data sets. The inversion results from synthetic data show the ability of characterizing anomalies of low- and high-velocity zones, and the inversion results from real data are generally consistent with SPT N-value, including the identification of a buried low-velocity layer.

Performances of the Partition of Unity Finite Element Method for the analysis of two-dimensional interior sound fields with absorbing materials

The paper deals with the numerical simulation of the acoustic field in two-dimensional cavities in which absorbing materials are present. Though Finite Element Method (FEM) could be employed for this purpose, the discretization level required for achieving reasonable accuracy renders the method impractical in the mid-frequency range. To alleviate this limitation, the Partition of Unity Finite Element Method (PUFEM) using plane wave functions has been shown to be very effective for solving short wave Helmholtz problems. In the present work, the method is extended to the computation of the pressure wave field within the absorbing media which is modeled as a bulk-reacting material characterized by a complex-valued and frequency dependent mean density and dynamic compressibility. Lagrange multipliers are used to enforce the transmission conditions at the air–material interface. Performances of the PUFEM are compared with a standard FEM in several examples of practical interests. It is shown that the technique is a good candidate for solving noise control problems at high frequency.

A SCENARIO STUDY FOR IMPROVING COST-EFFECTIVENESS IN ACOUSTIC TIME-REVERSAL SOURCE RELOCATION IN AN URBAN ENVIRONMENT

Through finite difference time domain (FDTD) numerical simulation, we have studied the possible observation settings to improve the cost effectiveness in time-reversal (TR) source relocation in a two-dimensional (2D) urban setting under a number of typical scenarios. All scenario studies were based on the FDTD computation of the acoustic wave field resulted from an impulse source, propagated through an artificial village composed of 15 buildings and a set of sources and receivers, a typical urban setting has been extensively analyzed in previous studies. The FDTD numerical modeling code can be executed on an off-the-shelf graphic processor unit (GPU) that increases the speed of the time-reversal calculations by a factor of 200. With this approach the computational results lead to some significant conclusions. In general, using only one non-line-of-sight (NLOS) single receiver is not enough to do a quality work to re-locate the source via time-reversal. This is particularly true when there are more than one path between the source and this receiver with similar wave energy travel time. However, when the single sensor is located in an acoustic channel, reverberation inside the waveguide may increase the effective aperture of the single receiver enough to give a good location. It is equivalent to say that the waveguide and the single receiver form a "virtual array". It appears that a sensor array with a minimum number of three receivers might be the most cost-effective way to carry out TR source relocation in an urban environment. The most optimal geometry of a sensor array with a minimum number of three receivers could be an equal side-length triangle. Simple analysis showed that by this setup it is possible to catch sound sources from almost all possible azimuths. Effective source relocation essentially depends on the geometry, relativity to the scatters, etc. of the sensing array. Generally, adding another single sensor relatively far away from the main array will not improve the results. It is practically useful and achievable to have a sensor array mounted on the outside of a single building, and in these cases successful source relocations were obtained. As stated by the fundamental TR theory, increasing the number of scatters, here, increasing the number of buildings will definitely be helpful to increase the effectiveness of TR source relocation.

Iterative method for the determination of a diffraction field in the interaction of a longitudinal shear wave with a system of cracks

UDC 539.3 We consider the problem of determination of a diffraction field that is formed due to the interaction of longitudinal shear waves with a system of cracks arbitrarily located in the infinite body. The initial problem is reduced to a system of singular integro-differential equations. We propose an iterative method of solution of this system, where the zero approximation is the solutions of integral equations for separate cracks. Modern methods of the mechanics of deformable bodies (method of potentials, method of discontinuous solutions, etc.) allow one to easily reduce the problem of determination of a wave field in the body with an arbitrary system of cracks to a system of integral equations. As a matter of fact, the interaction of longitudinal-shear waves (SH-waves) with an arbitrary system of cracks in the infinite elastic body was studied within the method of integral equations in [7, 11]. In [4], an analogous problem for radially positioned cracks was solved within the method of discontinuous solutions. These studies were continued in [3], where an elastic body under conditions of a plane deformation was considered. In these works, it is proposed to solve the systems of integrodifferential equations for the opening displacements of cracks and their derivatives in order to avoid hypersingular integrals. For such systems, an efficient procedure of numerical solution is available. To determine the wave fields in bodies with systems of cracks, the method of boundary elements is used more and more frequently in recent times. For example, this method was applied to the solution of two- and three-dimensional problems in [6, 8] and [9, 10], respectively. However, an analysis of these works shows that the reduction of initial boundary-value problems to boundary integral equations encounters serious difficulties, which consist of the necessity to numerically solve the systems of integral equations of high dimensions. The number of equations of such systems is proportional to the number of cracks, and the kernels of the equations themselves have singularities of high orders. For this reason, all authors restricted themselves to the case of several, most often two, cracks, while performing the numerical calculations of wave fields and stress intensity factors. Therefore, the problem of determination of a stressed state and diffraction fields in bodies with systems of cracks is not completely solved. In the present work, we propose an iterative approach to the determination of a diffraction field and the stress intensity factor in the case of the interaction of an SH-wave with a system of N cracks. As a result, there is no necessity to solve a system of N integro-differential equations, and, at each iteration, N separate equations corresponding to isolated cracks are solved. Statement of the Problem Let N through cracks be located in an elastic body in the state of plane deformation. In the plane Oxy , these cracks occupy intervals 2d k in length, whose middle points have coordinates (a k , b k ) , k = 1, 2, … , N

Framework for computing the spatial coherence effects of polycapillary x-ray optics

Despite the extensive use of polycapillary x-ray optics for focusing and collimating applications, there remains a significant need for characterization of the coherence properties of the output wavefield. In this work, we present the first quantitative computational method for calculation of the spatial coherence effects of polycapillary x-ray optical devices. This method employs the coherent mode decomposition of an extended x-ray source, geometric optical propagation of individual wavefield modes through a polycapillary device, output wavefield calculation by ray data resampling onto a uniform grid, and the calculation of spatial coherence properties by way of the spectral degree of coherence.

ГЕОДИНАМІКА

The program for calculations of seismic wave fields in media with free boundaries was developed. Was modeled an influence of a relief upon intensity in points on the uneven day surface.

THE ADJOINT METHOD APPLIED TO TIME-DISTANCE HELIOSEISMOLOGY

For a given misfit function, a specified optimality measure of a model, its gradient describes the manner in which one may alter properties of the system to march toward a stationary point. The adjoint method, arising from partial-differential-equation-constrained optimization, describes a means of extracting derivatives of a misfit function with respect to model parameters through finite computation. It relies on the accurate calculation of wavefields that are driven by two types of sources, namely, the average wave-excitation spectrum, resulting in the forward wavefield, and differences between predictions and observations, resulting in an adjoint wavefield. All sensitivity kernels relevant to a given measurement emerge directly from the evaluation of an interaction integral involving these wavefields. The technique facilitates computation of sensitivity kernels (Fréchet derivatives) relative to three-dimensional heterogeneous background models, thereby paving the way for nonlinear iterative inversions. An algorithm to perform such inversions using as many observations as desired is discussed.

SWAN Nested model의 현장 적용성 분석

최근 추세에 따라 수치실험은 높은 해상도와 정확도를 요구하고 있다. 일반적으로 파랑장 계산은 광역모형을 도입하고, 주된 해석 대상 영역의 지형적 복잡성을 반영할 수 있는 적절한 분해능의 상세격자를 통해 내부의 상세역 모형을 구성하게 된다. 여기서, 세역에서의 입사파는 광역에서의 결과를 토대로 추출하게 되는데, 이 때 결과의 연속성에서 문제를 가지게 된다. 대체로 이러한 문제점을 극복하기 위하여 가변격자체계와 스펙트럼 추출 등의 방법을 사용하게 되지만 관련 수치적 오류에 대한 분석과 검토가 충분하지 않은 것으로 보인다. 따라서, 본 연구에서는 최근 가장 일반적으로 사용되고 있는 SWAN모형의 Nested 모형에 대한 현장 적용성을 확인하고자 한다. 이를 위하여 각기 다른 조석환경하에 있는 두 개의 항에 대한 실험을 실시하고 이에 대한 결과를 비교 분석하였다. 분석결과, Nested 모형으로 구성된 광역과 세역 경계에서의 추출된 값이 거의 동일한 값을 가지며 각기 다른 조석환경하에서의 결과 또한 다르지 않다는 것으로 나타났다. 그러나, Nested 모형의 경계를 구성시 수심의 변화가 심한 곳을 피하고 비교적 수심이 깊은 해역까지를 포함시키는 것이 오차를 줄일 수 있음을 확인하였다. The recent trend for numerical experiment requires more higher resolution and accuracy. Generally, in the wave field calculation, it starts with a large region formulation first and follows by a separated detailed region formulation by more denser grids for the main interest area considering the geographical and bathymetrical variation. The wave fields resulted from the large region calculation is being introduced into the detail region calculation as the incident waves. In this process there exists a problem of continuity. In order to get over such problem, method of variable gridding system or spectrum sampling, etc., is being used. However, it seems not enough to examine and analyze the related numerical errors. Therefore, it is investigated in this study the field applicability of the most pervasive use of wave model, the nested SWAN model. For this purpose, we made model experiment for two coastal harbours with different tidal environment, and compared and analyzed the result. From the analysis, it was found that both the extracted values, near the boundaries of the large and detail region and the nested formulation of SWAN model, show almost the same and no different between those with different tidal environment conditions. However it is necessary for reducing the numerical errors to set the boundaries for the detailed region outside of the rapid bathymetric change and deeper region.

ON BASIC CONDITIONS FOR LONG-WAVE SIMULATIONS IN HARBORS BY THE BOUSSINESQ MODEL

For wave field calculations to estimate long-period ship motions in harbors, the so-called Boussinesq model seems to be widely used because of its versatility. However, it has not yet been clarified how large the sponge layer should be to simulate an open-sea condition, how long the computation should be continued to obtain a statistically stable spectrum, and what size of computational area should be considered. In this study, these basic conditions were investigated by conducting wave field simulations using the Boussinesq model. Useful information was obtained such as: (1) the thickness of the sponge layer can be much less than the thickness expected from one-dimensional calculations, (2) the data length to obtain a stable spectrum needs to be much longer than statistically expected 3 to 5 hours, and this is because the beat waves are produced by selectively amplified component waves near resonant frequencies in harbors.(3) the size of the computational area does not affect the results greatly and the distance 3 km between the end of the breakwater and the lateral sponge boundary is practically sufficient for the harbor mouth opening.

Tide simulation using the mild-slope equation with Coriolis force and bottom friction

Since the mild-slope equation was derived by Berkhoff (1972), the researchers considered various mechanism to simplify and improve the equation, which has been widely used for coastal wave field calculation. Recently, some scholars applied the mild-slope equation in simulating the tidal motion, which proves that the equation is capable to calculate the tide in actual terrain. But in their studies, they made a lot of simplifications, and did not consider the effects of Coriolis force and bottom friction on tidal wave. In this paper, the first-order linear mild-slope equations are deduced from Kirby mild-slope equation including wave and current interaction. Then, referring to the method of wave equations’ modification, the Coriolis force and bottom friction term are considered, and the effects of which have been performed with the radial sand ridges topography. Finally, the results show that the modified mild-slope equation can be used to simulate tidal motion, and the calculations agree well with the measurements, thus the applicability and validity of the mild-slope equation on tidal simulation are further proved.