Pseudospectral Time-domain Method Research Articles

Calibration of absorbing boundary layers for geoacoustic wave modeling in pseudo-spectral time-domain methods

Abstract. This paper develops a calibration methodology of the artificial absorbing techniques typically used by Fourier pseudo-spectral time-domain (PSTD) methods for geoacoustic wave simulations. Specifically, we consider the damped wave equation (DWE) that results from adding a dissipation term to the original wave equation, the sponge boundary layer (SBL) that applies an exponentially decaying factor directly to the wavefields, and finally, a classical split formulation of the perfectly matched layer (PML). These three techniques belong to the same family of absorbing boundary layers (ABLs), where outgoing waves and edge reflections are progressively damped across a grid zone of NABL consecutive layers. The ABLs used are compatible with a pure Fourier formulation of PSTD. We first characterize the three ABLs with respect to multiple sets of NABL and their respective absorption parameters for homogeneous media. Next, we validate our findings in heterogeneous media of increasing complexity, starting with a layered medium and finishing with the SEG/EAGE 3D salt model. Finally, we algorithmically compare the three PSTD–ABL methods in terms of memory demands and computational cost. An interesting result is that PML, despite outperforming the absorption of the other two ABLs for a given NABL value, requires approximately double the computational time. The parameter configurations reported in this paper can be readily used for PSTD simulations in other test cases, and the ABL calibration methodology can be applied to other wave propagation schemes.

Fourier Neural Operator Network for Fast Photoacoustic Wave Simulations

Simulation tools for photoacoustic wave propagation have played a key role in advancing photoacoustic imaging by providing quantitative and qualitative insights into parameters affecting image quality. Classical methods for numerically solving the photoacoustic wave equation rely on a fine discretization of space and can become computationally expensive for large computational grids. In this work, we applied Fourier Neural Operator (FNO) networks as a fast data-driven deep learning method for solving the 2D photoacoustic wave equation in a homogeneous medium. Comparisons between the FNO network and pseudo-spectral time domain approach were made for the forward and adjoint simulations. Results demonstrate that the FNO network generated comparable simulations with small errors and was orders of magnitude faster than the pseudo-spectral time domain methods (~26× faster on a 64 × 64 computational grid and ~15× faster on a 128 × 128 computational grid). Moreover, the FNO network was generalizable to the unseen out-of-domain test set with a root-mean-square error of 9.5 × 10−3 in Shepp–Logan, 1.5 × 10−2 in synthetic vasculature, 1.1 × 10−2 in tumor and 1.9 × 10−2 in Mason-M phantoms on a 64 × 64 computational grid and a root mean squared of 6.9 ± 5.5 × 10−3 in the AWA2 dataset on a 128 × 128 computational grid.

Numerical model of light propagation through Fabry-Perot etalons composed of interfaces with non-planar surface topography.

We present a model that calculates optical fields reflected and transmitted by a Fabry-Perot (FP) etalon composed of interfaces with non-planar surface topography. The model uses the Rayleigh-Rice theory, which predicts the fields reflected and transmitted by a single interface, to account for the non-planar surface topography of each interface. The Rayleigh-Rice theory is evaluated iteratively to account for all round trips that light can take within the FP etalon. The model predictions can then be used to compute Interferometer transfer function (ITF)s, by performing wavelength or angle resolved simulations enabling predictions of the bandwidth, peak transmissivity, and sensitivity of FP etalons. The model was validated against the Pseudospectral time-domain (PSTD) method, which resulted in good agreement. Since the model accuracy is expected to reduce as the Root mean square (RMS) of the topographic map increases, the error in the model's predictions was studied as a function of topographic map RMS. Finally, application of the model was exemplified by predicting the impact of roughness on ITFs and computing the changes in FP etalon transmissivity as cavity thickness is modulated by an ultrasonic wave.

Pseudo-Spectral Time-Domain Method for Subsurface Imaging with the Lunar Regolith Penetrating Radar

Recently and successfully, the Chang’E-5 (CE-5) lander was launched on a mission to bring 1.731 kg of lunar soil back to Earth. To investigate various compositions of lunar regolith, we apply the Lunar Regolith Penetrating Radar (LRPR) as the same scientific payload installed on the CE-5 lander. Based on the high-accuracy imaging technique, we achieve subsurface imaging to process LRPR-measured data collected from the lunar-like exploration tests in our laboratory. In this paper, we propose the pseudo-spectral time-domain (PSTD) method as the underlying code to implement the reverse-time migration (RTM) method and restore the uncertain subsurface area. With the significant advantage of lower spatial sampling density, the PSTD-RTM method not only saves major computational resources, but also rapidly confirms the object prediction in the effective imaging area. To further analyze the LRPR measured data, we employ the spectrum window to remove high- and low-frequency noise, and thus improve imaging visibility to some extent. The imaging results in this paper can prove the reliability and efficiency of the PSTD-RTM method for subsurface discoveries in planetary exploration.

Application of a Taylor series approximation to the Debye-Wolf integral in time-domain numerical electromagnetic simulations.

Finite-difference time-domain (FDTD) and pseudospectral time-domain (PSTD) methods are numerical electromagnetic simulation techniques that have been employed to perform rigorous simulations of broadband illuminations in several contexts. However, the computational cost of calculating the incident source fields introduced into the FDTD/PSTD grid can be considerable. In some cases, this can exceed the computational cost of what might be considered the principal part of the FDTD/PSTD algorithm, which calculates the spatial derivative of fields throughout the computational grid. In this paper, we analyze an existing method that has been used to approximate broadband illumination, which uses knowledge of the field only at a central frequency of the spectrum. We then present a new, to the best of our knowledge, approximation of the broadband illumination, which is more accurate, while remaining computationally tractable. Finally, we present some examples to verify the accuracy and efficiency of the new method and compare these results with the existing method.

Stability and accuracy of Engquist–Majda absorbing boundary condition for pseudo spectral time domain method

AbstractLagrange polynomials based Chebyshev Pseudo Spectral Time Domain (L‐CPSTD) method is an accurate time domain solver for electromagnetic problems. It utilizes the Lagrange interpolation polynomials to expand electromagnetic fields. This global interpolation, which uses all field values on all of the grid points to calculate the value at a single point, provides spectral accuracy. However, absorbing boundary conditions (ABCs) must be applied for open space problems. Engquist–Majda ABC is an important one due to its simplicity. Characteristic variables (CVs) can be used to implement the ABCs. In this article, for the first time, stability and accuracy of the Engquist–Majda ABC are proved by using the CVs in the L‐CPSTD solution of Maxwell's equations. The theoretical findings are verified by using the matrix eigenvalue method and the reflection coefficient in one and two dimensional open space examples. The numerical result is also validated by an analytical and an FDTD solution of a parallel‐plate waveguide problem. The efficiency of the proposed ABC is clearly shown.

Constant-Q wave propagation and compensation by pseudo-spectral time-domain methods

Seismic attenuation is usually described by a constant-Q (CQ) model that assumes the seismic quality factor (Q) is independent of the frequency. To simulate the attenuation behaviors of seismic waves, we develop a pseudo-spectral time-domain (PSTD) method to solve a CQ viscoacoustic wave equation. This method is nearly fourth-order accurate in time. Compared to the conventional temporal second-order PSTD method, the new PSTD scheme is verified to be more efficient. In some applications such as Q-compensated reverse-time migration (Q-RTM) and time-reversal imaging, one requires to simulate an anti-attenuation process. To realize this purpose, we switch our viscoacoustic PSTD modeling scheme into an amplitude-compensated PSTD modeling scheme by flipping the signs of the operators that dominate the amplitude loss. To control the numerical instability caused by high-frequency overcompensation in the amplitude-boosted modeling, we integrate a time-variant filter to the PSTD modeling scheme to suppress the high-frequency noise. Wavefield simulation examples in homogeneous media verify the temporal accuracy of our nearly fourth-order PSTD modeling scheme. A Q-RTM test of synthetic data is also presented to demonstrate the robustness of our amplitude-compensated PSTD modeling scheme.

A Generalized Split-Step Angular Spectrum Method for Efficient Simulation of Wave Propagation in Heterogeneous Media.

Angular spectrum (AS) methods enable efficient calculation of wave propagation from one plane to another inside homogeneous media. For wave propagation in heterogeneous media such as biological tissues, AS methods cannot be applied directly. Split-stepping techniques decompose the heterogeneous domain into homogeneous and perturbation parts, and provide a solution for forward wave propagation by propagating the incident wave in both frequency-space and frequency-wavenumber domains. Recently, a split-step hybrid angular spectrum (HAS) method was proposed for plane wave propagation of focused ultrasound beams. In this study, we extend these methods to enable simulation of acoustic pressure field for an arbitrary source distribution, by decomposing the source and reflection spectra into orthogonal propagation direction components, propagating each component separately, and summing all components to get the total field. We show that our method can efficiently simulate the pressure field of arbitrary sources in heterogeneous media. The accuracy of the method was analyzed comparing the resultant pressure field with pseudospectral time domain (PSTD) solution for breast tomography and hemispherical transcranial-focused ultrasound simulation models. Eighty times acceleration was achieved for a 3-D breast simulation model compared to PSTD solution with 0.005 normalized root mean-squared difference (NRMSD) between two solutions. For the hemispherical phased array, aberrations due to skull were accurately calculated in a single simulation run as evidenced by the resultant-focused ultrasound beam simulations, which had 0.001 NRMSD with 40 times acceleration factor compared to the PSTD method.

A Simple Solver for the Two-Fluid Plasma Model Based on PseudoSpectral Time-Domain Algorithm

We present a solver of 3D two-fluid plasma model for the simulation of short-pulse laser interactions with plasma. This solver resolves the equations of the two-fluid plasma model with ideal gas closure. We also include the Bhatnagar-Gross-Krook collision model. Our solver is based on PseudoSpectral Time-Domain (PSTD) method to solve Maxwell's curl equations. We use a Strang splitting to integrate Euler equations with source term: while Euler equations are solved with a composite scheme mixing Lax-Friedrichs and Lax-Wendroff schemes, the source term is integrated with a fourth-order Runge-Kutta scheme. This two-fluid plasma model solver is simple to implement because it only relies on finite difference schemes and Fast Fourier Transforms. It does not require spatially staggered grids. The solver was tested against several well-known problems of plasma physics. Numerical simulations gave results in excellent agreement with analytical solutions or with previous results from the literature.

Reverse-Time Migration by Combining Laplacian Filtering With Wavefield Decomposition for High-Resolution Subsurface Imaging

Reverse-time migration (RTM) is considered to be one of the most accurate migration methods, especially for imaging geologically complex structures. However, the conventional RTM using cross-correlation imaging condition is subject to strong migration artifacts. Recent researches on electromagnetic wavefield decomposition exhibit that the method can greatly suppresses the internal reflection noise, however, it is still subject to residual noise. To further eliminate the low-frequency noise to obtain high-resolution subsurface image, a novel electromagnetic imaging scheme by combining a Laplacian filtering with wavefield decomposition is presented, and the pseudospectral time-domain (PSTD) method is introduced to solve the time-dependent partial differential equations efficiently. In this work, we apply the Laplacian filtering on each wavefield snapshot in space domain to sharpen the field distribution. The filtered wavefields are further decomposed into downgoing and upgoing components for final imaging. Numerical experiments show that the proposed method balances the amplitudes of reflectors, and exhibits higher resolution compared to the conventional RTM.

Full-Matrix Phase Shift Migration Method for Transcranial Ultrasonic Imaging.

A spectrum-domain method, called full-matrix phase shift migration (FM-PSM), is presented for transcranial ultrasound phase correction and imaging with ideal synthetic aperture focusing technology. The simulated data obtained using the pseudospectral time-domain method are used to evaluate the feasibility of the method. The experimental data measured from a 3-D printed skull phantom are used to evaluate the algorithm performance in terms of resolution, contrast-to-noise ratio (CNR), and eccentricity comparing with the classical ray-tracing delay and sum (DAS) method. In wire imaging experiment, FM-PSM has a lateral resolution of 0.22 mm and ray-tracing DAS has a lateral resolution of 0.24 mm measured at -6-dB drop using a transducer with a center frequency of 6.25 MHz. In cylinder imaging experiment, FM-PSM has a CNR of 2.14 and ray-tracing DAS has a CNR of 1.82, which illustrates about 17% improvement. For a J -element array and an output image with pixels M ×N (lateral × axial), the computational cost of the DAS is of O(J ×M2×N2) ; on the contrary, it reduces to O(J ×M ×N2) with the proposed FM-PSM. The results suggest that FM-PSM is an efficiency method for transcranial ultrasonic imaging.

An efficient and accurate parallel hybrid acoustic signal correction method for transcranial ultrasound

Phased-control focused ultrasound transducers provide a new and noninvasive treatment method for brain disease. However, improving the accuracy of phase correction and reducing the calculation time during treatment have always been contradictory constraints. In this paper, a hybrid acoustic signal correction (HASC) method combined with k-Wave stage and holography stage was introduced for phase correction and simulation of transcranial focused ultrasound. The k-Wave stage is mainly used to calculate the sound field in a heterogeneous medium (skull), which divides the sound field calculation process into paths that can be calculated in parallel, and the transcranial correction phase can also be obtained during the calculation. The holography stage is sufficient to simulate the acoustic field in the homogenous intracranial medium after ultrasound transmitting through the skull. The agreement of the k-space corrected pseudospectral time domain method and HASC method was assessed by statistical methods: linear regression between the two methods provided a slope of 0.9735, intercept of 0.0078, and R2 of 0.9982. The Bland–Altman method provided a bias of 0.0015 and 95% limits of agreement 0.065 apart. We demonstrated that the difference in sound intensity at the focal point corrected by HASC and time reversal phase correction method was 0.2% and 0.5% in the results of simulation and experiment, respectively. Not only that, the phase calculation time by the HASC phase correction method can be reduced to 11 min on a multi GPU array, which has clinical potential for ultrasound treatment of brain therapy.

Pseudospectral Time-Domain (PSTD) Methods for the Wave Equation: Realizing Boundary Conditions with Discrete Sine and Cosine Transforms

Pseudospectral time domain (PSTD) methods are widely used in many branches of acoustics for the numerical solution of the wave equation, including biomedical ultrasound and seismology. The use of the Fourier collocation spectral method in particular has many computational advantages. However, the use of a discrete Fourier basis is also inherently restricted to solving problems with periodic boundary conditions. Here, a family of spectral collocation methods based on the use of a sine or cosine basis is described. These retain the computational advantages of the Fourier collocation method but instead allow homogeneous Dirichlet (sound-soft) and Neumann (sound-hard) boundary conditions to be imposed. The basis function weights are computed numerically using the discrete sine and cosine transforms, which can be implemented using [Formula: see text] operations analogous to the fast Fourier transform. Practical details of how to implement spectral methods using discrete sine and cosine transforms are provided. The technique is then illustrated through the solution of the wave equation in a rectangular domain subject to different combinations of boundary conditions. The extension to boundaries with arbitrary real reflection coefficients or boundaries that are nonreflecting is also demonstrated using the weighted summation of the solutions with Dirichlet and Neumann boundary conditions.

Nonlinear ultrasound simulation in an axisymmetric coordinate system using a k-space pseudospectral method.

A full-wave model for nonlinear ultrasound propagation through a heterogeneous and absorbing medium in an axisymmetric coordinate system is developed. The model equations are solved using a nonstandard or k-space pseudospectral time domain method. Spatial gradients in the axial direction are calculated using the Fourier collocation spectral method, and spatial gradients in the radial direction are calculated using discrete trigonometric transforms. Time integration is performed using a k-space corrected finite difference scheme. This scheme is exact for plane waves propagating linearly in the axial direction in a homogeneous and lossless medium and significantly reduces numerical dispersion in the more general case. The implementation of the model is described, and performance benchmarks are given for a range of grid sizes. The model is validated by comparison with several analytical solutions. This includes one-dimensional absorption and nonlinearity, the pressure field generated by plane-piston and bowl transducers, and the scattering of a plane wave by a sphere. The general utility of the model is then demonstrated by simulating nonlinear transcranial ultrasound using a simplified head model.

Performance analysis of parallelized PSTD-FDTD method for large-scale electromagnetic simulation

We developed a parallel solver for large-scale electromagnetic simulation by hybridizing the Pseudo-Spectral Time Domain (PSTD) and Finite Difference Time Domain (FDTD) methods. Parallelization is implemented with a Message Passing Interface (MPI) and Open Multi-Processing (OpenMP) and applied to a single spatial dimension in the full three-dimensional simulation. The stability and accuracy of the hybridized PSTD-FDTD (HPF) method are analyzed and compared with the FDTD method, in particular by applying to a known exact solution of the Mie scattering theory. We also present a roofline model analysis to explain the high computing speed of the HPF method despite its heavy computational loads. The HPF method features high stability, low parallelization cost and fast computing speed on the PC cluster, especially when the sampling rate is low.

Stabilized DG-PSTD Method With Nonconformal Meshes for Electromagnetic Waves

We present a node-based discontinuous Galerkin (DG) pseudospectral time domain (PSTD) algorithm, with adaptive nonconformal unstructured meshes, for 3-D large-scale Maxwell’s equations. This algorithm is a combination of a new DG algorithm and a PSTD method, where spectral accuracy is achieved via the PSTD algorithm, while the DG serves as a stable coupling for multiple domains with unstructured hexahedra. Time marching is efficient because the mass matrix in the DG-PSTD algorithm is exactly diagonal. The scheme is low-storage and scalable because the stiffness matrix is localized into a small shared matrix. Furthermore, arbitrary nonconformal meshes can be adaptively realized, increasing the flexibility of complex media modeling. Our numerical results corroborate the long-time stability, high efficiency, and high-order accuracy of the proposed solver. Finally, an adaptive application of 5G electromagnetic signal propagation demonstrates the efficiency and capability of the proposed high-order solver.

Scattering matrices of mineral dust aerosols: a refinement of the refractive index impact

Abstract. Mineral dust, as one of the most important aerosols, plays a crucial role in the atmosphere by directly interacting with radiation, while there are significant uncertainties in determining dust optical properties to quantify radiative effects and to retrieve their properties. Laboratory and in situ measurements of the refractive indices (RIs) of dust differ, and different RIs have been applied in numerical studies used for model developments, aerosol retrievals, and radiative forcing simulations. This study reveals the importance of the dust RI for the development of a model of dust optical properties. The Koch-fractal polyhedron is used as the modeled geometry, and the pseudospectral time domain method and improved geometric-optics method are combined for optical property simulations over the complete size range. We find that the scattering matrix elements of different kinds of dust particles are reasonably reproduced by choosing appropriate RIs, even when using a fixed particle geometry. The uncertainty of the RI would greatly affect the determination of the geometric model, as a change in the RI, even in the widely accepted RI range, strongly affects the shape parameters used to reproduce the measured dust scattering matrix elements. A further comparison shows that the RI influences the scattering matrix elements in a different way than geometric factors, and, more specifically, the P11, P12, and P22 elements seem more sensitive to the RI of dust. In summary, more efforts should be devoted to account for the uncertainties on the dust RI in modeling its optical properties, and the development of corresponding optical models can potentially be simplified by considering only variations over different RIs. Considerably more research, especially from direct measurements, should be carried out to better constrain the uncertainties related to the dust aerosol RIs.

Long-Time Instability of Pseudospectral Time-Domain Method in Curvilinear Coordinates

We found that the pseudospectral time-domain (PSTD) formulation discussed in the literature exhibits long-time instability when applied to problems with curvilinear geometry. To overcome this problem, we have developed a provably stable PSTD formulation in curvilinear coordinates. In this article, we discuss the new PSTD formulation and present a rigorous proof of stability. Numerical results are given to validate the new formulation. The new formulation requires only a simple change in the form of the spatial-derivative operators in curvilinear coordinates and so can easily be implemented in the existing PSTD codes.

Auralization of a Car Pass-By Using Impulse Responses Computed with a Wave-Based Method

Auralization of traffic noise can be of great value for city planners and in the communication to citizens. In inner city environments, auralization of traffic noise based on predictions might require a wave-based method for the low frequencies. In this paper, a method to auralize a car pass-by using binaural impulse responses (BIRs) computed with the wave-based pseudospectral time-domain method (PSTD) in two dimensions is presented. A dry synthesized car signal is convolved with the BIRs for different locations in the street, and cross-fade windows are used in order to create a smooth transition between the source positions. Auralizations were performed for the simplified scenarios where buildings are absent, and for an environment where a long flat wall is located behind the car, opposite from the receiver. A same-different listening test was carried out in order to investigate whether increasing the angular spacing between the discrete source positions affects the perception of the car pass-by auralizations. Signal detection theory (SDT) was used for the design and analysis of the listening test. Results showed that differences exist, although they are difficult to notice. On average, 52.3% of the subjects found it almost impossible to spot any difference between auralizations with larger angular spacing (up to 10°) compared to the reference auralization (2° angular spacing).

Fractional Laplacian Viscoacoustic Wave Equation Low-Rank Temporal Extrapolation

Description of seismic wave attenuation is a hot topic in the geophysical area, and it is the basis of the attenuation-compensated seismic imaging technique, which aims to retrieve a high-resolution subsurface image for geological structure analysis and hydrocarbon reservoir prediction. The numerical simulation of viscoacoustic wave equation is an effective way to observe the seismic attenuation in lossy media. The existing study confirms that the seismic-quality-factor (Q) that is used to represent the strength of the viscous behavior of the earth is nearly independent of frequency, which is referred to as the constant-Q (CQ) model. The mathematical concept of fractional Laplacian is recently introduced to the geophysical area to form a compact CQ wave equation to describe seismic wave propagation. However, numerically solving the fractional Laplacian CQ wave equation by the traditional pseudospectral time-domain (PSTD) method suffers from a strict stability condition and great numerical dispersion due to a low-order temporal finite-difference (FD) approximation. To improve the temporal extrapolation accuracy, we derive the analytical wavenumber (k)-space domain propagators underlying the fractional Laplacian wave equation. We regard the k-space operators as mixed-domain matrices in the case of heterogeneous media and adopt a low-rank decomposition to approximate the matrices. With the low-rank approximation, an efficient time-marching formula is obtained for wavefield temporal extrapolation. We formulate the time-marching formula into a first-order equation system in terms of pressure and particle-velocity to welcome the perfectly matched layer (PML) absorbing boundary condition to eliminate the wraparound effects caused by the Fourier transform. A spatial-variable density is also incorporated to simulate more realistic amplitude variation. The numerical examples are carried out to verify the accuracy and stability of the viscoacoustic low-rank extrapolation.