Finite-difference Stencil Research Articles

Computationally efficient three-dimensional acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver

The computational burden of frequency-domain full-waveform inversion (FWI) of wide-aperture fixed-spread data is conventionally reduced by limiting the inversion to a few discrete frequencies. In this framework, frequency-domain seismic modeling is performed efficiently for multiple sources by solving the linear system resulting from the discretization of the time-harmonic wave equation with the massively parallel sparse direct solver. Frequency-domain seismic modeling based on the sparse direct solver (DSFDM) requires specific design finite-difference stencils of compact support to minimize the computational cost of the lower-upper decomposition of the impedance matrix in terms of memory demand and floating-point operations. A straightforward adaptation of such finite-difference stencil, originally developed for the (isotropic) acoustic-wave equation, is proposed to introduce vertical transverse isotropy (VTI) in the modeling without any extra computational cost. The method relies on a fourth-order wave equation, which is decomposed as the sum of a second-order elliptic-wave equation plus an anellipticity correction term. The stiffness matrix of the elliptic-wave equation is easily built from the isotropic stiffness matrix by multiplying its coefficients with factors that depend on Thomsen’s parameters, whereas the anelliptic term is discretized with a parsimonious second-order staggered-grid stencil. Validation of DSFDM against finite-difference time-domain modeling performed in various synthetic models shows that a discretization rule of four grid points per minimum wavelength provides accurate DSFDM solutions. Moreover, comparison between real data from the Valhall field and DSFDM solutions computed in a smooth VTI subsurface model supports that the method can be used as a fast and accurate modeling engine to perform multiparameter VTI FWI of fixed-spread data in the viscoacoustic approximation.

Reverse time migration from topography

Migration of seismic data from topography using methods based on finite-difference (FD) approximation to acoustic wave propagation commonly suffers from a number of imaging drawbacks due to the difficulty of applying FD stencils to irregular computational meshes. Altering the computational geometry from Cartesian to a topographic coordinate system conformal to the data acquisition surface can circumvent many of these issues. The coordinate transformation approach allows for acoustic wave propagation and the crosscorrelation and inverse-scattering imaging conditions to be posed and computed directly in topographic coordinates. Resulting reverse time migration (RTM) images may then be interpolated back to the Cartesian domain using the known inverse mapping. Orthogonal 2D topographic coordinates can be developed using known conformal mapping transforms and serve as the computational mesh for performing migration from topography. Impulse response tests demonstrate the accuracy of the 2D generalized acoustic wave propagation. RTM imaging examples show the efficacy of performing migration from topography directly from the data acquisition surface on topographic meshes and the ability to image complex near-surface structure even in the presence of strong lateral velocity variation.

Modeling of pseudoacoustic P-waves in orthorhombic media with a low-rank approximation

Wavefield extrapolation in pseudoacoustic orthorhombic anisotropic media suffers from wave-mode coupling and stability limitations in the parameter range. We use the dispersion relation for scalar wave propagation in pseudoacoustic orthorhombic media to model acoustic wavefields. The wavenumber-domain application of the Laplacian operator allows us to propagate the P-waves exclusively, without imposing any conditions on the parameter range of stability. It also allows us to avoid dispersion artifacts commonly associated with evaluating the Laplacian operator in space domain using practical finite-difference stencils. To handle the corresponding space-wavenumber mixed-domain operator, we apply the low-rank approximation approach. Considering the number of parameters necessary to describe orthorhombic anisotropy, the low-rank approach yields space-wavenumber decomposition of the extrapolator operator that is dependent on space location regardless of the parameters, a feature necessary for orthorhombic anisotropy. Numerical experiments that the proposed wavefield extrapolator is accurate and practically free of dispersion. Furthermore, there is no coupling of qSv and qP waves because we use the analytical dispersion solution corresponding to the P-wave.

Validation of adaptive rectangular decomposition for three-dimensional wave-based acoustic simulation in architectural models

Computer-based simulation is an increasingly popular way to predict the acoustics of real-world architectural designs. Most commercial acoustic simulation tools are based on geometric techniques and cannot accurately model low-frequency diffraction and other wave phenomena. Numerical wave simulation techniques can model these effects, but are less commonly used, since they are compute- and memory-intensive, and cannot scale to large spaces. Moreover, it is challenging to ensure that numerical methods do not suffer from high dispersion errors. Recent techniques have begun to overcome these limitations. One such method is adaptive rectangular decomposition (ARD), which combines analytical solutions to the wave equation in rectangular subdomains with a finite difference stencil for interface handling between subdomains, resulting in high-performance wave simulation with low dispersion error. ARD, along with high-performance ray tracing, are available as part of Impulsonic’s IPL SDK, a software development kit that allows custom acoustic simulation tools to be easily built with state-of-the-art simulation technology. In this paper, we evaluate the performance and accuracy of the IPL SDK and ARD, by analyzing simulation results and comparing them against measurements obtained for real-world architectural designs.

Efficient parallel implementation of incompressible pipe flow algorithm based on SIMPLE

SummaryParallel semi‐implicit method for pressure‐linked equations(SIMPLE) algorithm is used to solve the 3‐D incompressible pipe flow problem. In this paper, we proposed a novel parallel SIMPLE algorithm that uses the alternate tiling technique. Firstly, a parallel SIMPLE algorithm based on domain decomposition method was established, and the implementation of domain partition and data exchange was presented. Then, we presented serial finite difference stencil algorithm based on alternate tiling. Furthermore, an iteration space parallel two‐way finite difference stencil algorithm based on alternate tiling was proposed, introducing the sequence of iterative space tiles as the sequence of execution and using time skewing technique to partition the iteration space, thus to improve the data locality of algorithm. The cache misses and the cost of communication and synchronization are reduced by reordering the tiles of iteration space. Finally, the effectiveness of the two parallel SIMPLE algorithms were compared. The results showed that the parallel SIMPLE algorithm that uses the two‐way finite difference stencil algorithm based on alternate tiling has good data locality, performance, and scalability in the Deepcomp7000 cluster computing environment. Copyright © 2013 John Wiley & Sons, Ltd.

A novel concept for the design of immersed interface methods

The objective of this paper is to present a novel, robust, high-order accurate Immersed Interface Method (IIM) for advection–diffusion type equations. In contrast to other immersed methods that were designed for consistency and accuracy with a posteriori check of the numerical stability, we combine local Taylor-series expansion at irregular grid points with a local stability constraint as part of the design process. Stability investigations of the IIM are employed to demonstrate that the local stability constraint is sufficient for obtaining a globally stable method, as long as the Neumann number is less than its limiting value. One of the key aspects of this IIM is that the irregular finite-difference stencils can be isolated from the rest of the computational domain. To validate our novel immersed interface approach, two-dimensional and three-dimensional test cases for model equations are presented. In addition, this method is applied to the incompressible Navier–Stokes equations to conduct stability investigations of a boundary layer flow over a rough surface, and for investigations of pulsatile stenotic flows. Stability investigations of wall bounded flows are challenging for immersed methods, because the near wall accuracy is important for correctly capturing the characteristics of the hydrodynamic instability mechanisms, in particular regarding the wave relation between the wave velocity components close to the wall.

Recursive evaluation of space-time lattice Green's functions

Up to a multiplicative constant, the lattice Green's function (LGF) as defined in condensed matter physics and lattice statistical mechanics is equivalent to the Z-domain counterpart of the finite-difference time-domain Green's function (GF) on a lattice. Expansion of a well-known integral representation for the LGF on a ν-dimensional hyper-cubic lattice in powers of Z−1 and application of the Chu–Vandermonde identity results in ν − 1 nested finite-sum representations for discrete space-time GFs. Due to severe numerical cancellations, these nested finite sums are of little practical use. For ν = 2, the finite sum may be evaluated in closed form in terms of a generalized hypergeometric function. For special lattice points, that representation simplifies considerably, while on the other hand the finite-difference stencil may be used to derive single-lattice-point second-order recurrence schemes for generating 2D discrete space-time GF time sequences on the fly. For arbitrary symbolic lattice points, Zeilberger's algorithm produces a third-order recurrence operator with polynomial coefficients of the sixth degree. The corresponding recurrence scheme constitutes the most efficient numerical method for the majority of lattice points, in spite of the fact that for explicit numeric lattice points the associated third-order recurrence operator is not the minimum recurrence operator. As regards the asymptotic bounds for the possible solutions to the recurrence scheme, Perron's theorem precludes factorial or exponential growth. Along horizontal lattices directions, rapid initial growth does occur, but poses no problems in augmented dynamic-range fixed precision arithmetic. By analysing long-distance wave propagation along a horizontal lattice direction, we have concluded that the chirp-up oscillations of the discrete space-time GF are the root cause of grid dispersion anisotropy. With each factor of ten increase in the lattice distance, one would have to roughly double the pulse width of the source signature to keep pulse distortion at bay. The GF time sequences can also be used for an efficient computation of discrete space-frequency LGFs, especially if one employs Aitken's δ2 process for the acceleration of the convergence of the consecutive partial sums.

Time–space domain dispersion-relation-based finite-difference method with arbitrary even-order accuracy for the 2D acoustic wave equation

High-order finite-difference (FD) methods have been widely used for numerical solution of acoustic wave equations. It has been reported that the modeling accuracy is of 2nd-order when the conventional (2M)th-order space domain FD and the 2nd-order time domain FD stencils are directly used to solve the acoustic wave equation. Under the same discretization, the present version of time–space domain dispersion-relation-based FD method can improve the accuracy from 2nd-order to (2M)th-order along eight directions for the 2D acoustic wave equation. To increase the accuracy further, we propose a new FD stencil for 2D acoustic wave equation modeling. This new time–space domain dispersion-relation-based FD stencil can reach the same arbitrary even-order accuracy along all directions, and is more accurate and more stable than the conventional one for the same M. Dispersion analysis and modeling examples demonstrate its advantages.

Sheath governing equations in computational weakly-ionized plasmadynamics

To date, fluid models of plasma sheaths have consisted of the coupling of the electric field potential equation obtained through Gauss’s law to the charged species conservation equations obtained through the drift–diffusion approximation. When discretized using finite-difference stencils, such a set of equations has been observed to be particularly stiff and to often require more than hundreds of thousands of iterations to reach convergence. A new approach at solving sheaths using a fluid model is here presented that reduces significantly the number of iterations to reach convergence while not sacrificing on the accuracy of the converged solution. The method proposed herein consists of rewriting the sheath governing equations such that the electric field is obtained from Ohm’s law rather than from Gauss’s law. To ensure that Gauss’s law is satisfied, some source terms are added to the ion conservation equation. Several time-accurate and steady-state cases of dielectric sheaths, anode sheaths, and cathode sheaths (including glow and dark discharges) are considered. The proposed method is seen to yield the same converged solution as the conventional approach while exhibiting a reduction in computational effort varying between one-hundred-fold and ten-thousand-fold whenever the plasma includes both quasi-neutral regions and non-neutral sheaths.

An iterated, multipoint differential transform method for numerically evolving partial differential equation initial-value problems

Traditional numerical techniques for solving time-dependent partial differential equation (PDE) initial-value problems (IVPs) store a truncated representation of the function values and a certain number of their time derivatives at each time step. Although redundant in the dx → 0 limit, what if spatial derivatives were also stored? This paper presents an iterated, multipoint differential transform method (IMDTM) for numerically evolving PDE IVPs. Using this scheme, we demonstrate that stored spatial derivatives can be propagated in an efficient and self-consistent manner and can effectively contribute to the evolution procedure in a way that can confer several advantages, including aiding in solution verification. Lastly, in order to efficiently implement the IMDTM scheme, a generalized finite-difference stencil formula is derived that can take advantage of multiple higher-order spatial derivatives when computing even-higher-order derivatives. As demonstrated here, the performance of these techniques compares favorably to other explicit evolution schemes in terms of speed, memory footprint and accuracy.

3D scalar-wave absorbing boundary conditions with optimal coefficients in the frequency-space domain

To improve the computational efficiency for the solution of the 3D Helmholtz equation in the frequency-space domain, high-order compact forms of finite differences are preferred. We applied a pointwise Padé approximation to develop a 3D 27-point fourth-order compact finite-difference (FD) stencil in the grid interior, with a space-differentiated source term, for the scalar-wave equation; this has similar high-accuracy (4–5 grid points per the shortest wavelength) to another 27-point fourth-order FD stencil using a parsimonious mixed-grid and staggered-grid combination, but is much simpler. For absorbing boundary conditions (ABCs), a damping zone is expensive, and a perfectly matched layer can not be straightforwardly introduced into the compact FD form for the second-order wave equation. Thus, we developed 3D one-way wave equation (OWWE) ABCs with adjustable coefficients. They have different angle approximations and FD forms for the six faces, twelve edges, and eight corners in 3D models to fit with the interior compact FD form. By adjusting the coefficients to the optimum, the OWWE ABCs have wider-angle absorbing ability than those without optimal coefficients. Finally, all the interior and boundary FD forms were combined into a sparse complex-valued impedance matrix of the frequency-space modeling equation, and solved for each frequency. Because the storage of the sparse impedance matrix was determined by the 3D discrete grid size, the OWWE ABCs with only one outer layer needed the minimum grid size compared with other ABCs, thus were the most efficient for the solution of the impedance matrix. The modeling algorithm was performed on multicore processors using a MPI parallel direct solver. Numerical tests on homogeneous and heterogeneous models gave satisfactory absorbing effects.

Truncation Errors, Exact And Heuristic Stability Analysis Of Two-Relaxation-Times Lattice Boltzmann Schemes For Anisotropic Advection-Diffusion Equation

AbstractThis paper establishes relations between the stability and the high-order truncated corrections for modeling of the mass conservation equation with the two-relaxation-times (TRT) collision operator. First we propose a simple method to derive the truncation errors from the exact, central-difference type, recurrence equations of the TRT scheme. They also supply its equivalent three-time-level discretization form. Two different relationships of the two relaxation rates nullify the third (advection) and fourth (pure diffusion) truncation errors, for any linear equilibrium and any velocity set. However, the two relaxation times alone cannot remove the leading-order advection-d if fusion error, because of the intrinsic fourth-order numerical diffusion. The truncation analysis is carefully verified for the evolution of concentration waves with the anisotropic diffusion tensors. The anisotropic equilibrium functions are presented in a simple but general form, suitable for the minimal velocity sets and the d2Q9, d3Q13, d3Q15 and d3Q19 velocity sets. All anisotropic schemes are complemented by their exact necessary von Neumann stability conditions and equivalent finite-difference stencils. The sufficient stability conditions are proposed for the most stable (OTRT) family, which enables modeling at any Peclet numbers with the same velocity amplitude. The heuristic stability analysis of the fourth-order truncated corrections extends the optimal stability to larger relationships of the two relaxation rates, in agreement with the exact (one-dimensional) and numerical (multi-dimensional) stability analysis. A special attention is put on the choice of the equilibrium weights. By combining accuracy and stability predictions, several strategies for selecting the relaxation and free-tunable equilibrium parameters are suggested and applied to the evolution of the Gaussian hill.

3D acoustic wave modelling with time-space domain dispersion-relation-based finite-difference schemes and hybrid absorbing boundary conditions

Most conventional finite-difference methods adopt second-order temporal and (2M)th-order spatial finite-difference stencils to solve the 3D acoustic wave equation. When spatial finite-difference stencils devised from the time-space domain dispersion relation are used to replace these conventional spatial finite-difference stencils devised from the space domain dispersion relation, the accuracy of modelling can be increased from second-order along any directions to (2M)th-order along 48 directions. In addition, the conventional high-order spatial finite-difference modelling accuracy can be improved by using a truncated finite-difference scheme. In this paper, we combine the time-space domain dispersion-relation-based finite difference scheme and the truncated finite-difference scheme to obtain optimised spatial finite-difference coefficients and thus to significantly improve the modelling accuracy without increasing computational cost, compared with the conventional space domain dispersion-relation-based finite difference scheme. We developed absorbing boundary conditions for the 3D acoustic wave equation, based on predicting wavefield values in a transition area by weighing wavefield values from wave equations and one-way wave equations.Dispersion analyses demonstrate that high-order spatial finite-difference stencils have greater accuracy than low-order spatial finite-difference stencils for high frequency components of wavefields, and spatial finite-difference stencils devised in the time-space domain have greater precision than those devised in the space domain under the same discretisation. The modelling accuracy can be improved further by using the truncated spatial finite-difference stencils. Stability analyses show that spatial finite-difference stencils devised in the time-space domain have better stability condition. Numerical modelling experiments for homogeneous, horizontally layered and Society of Exploration Geophysicists/European Association of Geoscientists and Engineers salt models demonstrate that this modelling scheme has greater accuracy than a conventional scheme and has better absorbing effects than Clayton-Engquist absorbing boundary conditions.

An axis treatment for flow equations in cylindrical coordinates based on parity conditions

A novel axis treatment using parity conditions is presented for flow equations in cylindrical coordinates that are represented in azimuthal Fourier modes. The correct parity states of scalars and the velocity vector are derived such that symmetry conditions for each Fourier mode of the respective variable can be determined. These symmetries are then used to construct finite-difference and filter stencils at and near the axis, and an interpolation scheme for the computation of terms premultiplied by 1/ r. A grid convergence study demonstrates that the axis treatment retains the formal accuracy of the spatial discretization scheme employed. Two further test cases are considered for evaluation of the axis treatment, the propagation of an acoustic pulse and direct numerical simulation of a fully turbulent supersonic axisymmetric wake. The results demonstrate the applicability of the axis treatment for non-axisymmetric flows.

A Fast and Robust Poisson–Boltzmann Solver Based on Adaptive Cartesian Grids

An adaptive Cartesian grid (ACG) concept is presented for the fast and robust numerical solution of the 3D Poisson-Boltzmann Equation (PBE) governing the electrostatic interactions of large-scale biomolecules and highly charged multi-biomolecular assemblies such as ribosomes and viruses. The ACG offers numerous advantages over competing grid topologies such as regular 3D lattices and unstructured grids. For very large biological molecules and multi-biomolecule assemblies, the total number of grid-points is several orders of magnitude less than that required in a conventional lattice grid used in the current PBE solvers thus allowing the end user to obtain accurate and stable nonlinear PBE solutions on a desktop computer. Compared to tetrahedral-based unstructured grids, ACG offers a simpler hierarchical grid structure, which is naturally suited to multigrid, relieves indirect addressing requirements and uses fewer neighboring nodes in the finite difference stencils. Construction of the ACG and determination of the dielectric/ionic maps are straightforward, fast and require minimal user intervention. Charge singularities are eliminated by reformulating the problem to produce the reaction field potential in the molecular interior and the total electrostatic potential in the exterior ionic solvent region. This approach minimizes grid-dependency and alleviates the need for fine grid spacing near atomic charge sites. The technical portion of this paper contains three parts. First, the ACG and its construction for general biomolecular geometries are described. Next, a discrete approximation to the PBE upon this mesh is derived. Finally, the overall solution procedure and multigrid implementation are summarized. Results obtained with the ACG-based PBE solver are presented for: (i) a low dielectric spherical cavity, containing interior point charges, embedded in a high dielectric ionic solvent - analytical solutions are available for this case, thus allowing rigorous assessment of the solution accuracy; (ii) a pair of low dielectric charged spheres embedded in a ionic solvent to compute electrostatic interaction free energies as a function of the distance between sphere centers; (iii) surface potentials of proteins, nucleic acids and their larger-scale assemblies such as ribosomes; and (iv) electrostatic solvation free energies and their salt sensitivities - obtained with both linear and nonlinear Poisson-Boltzmann equation - for a large set of proteins. These latter results along with timings can serve as benchmarks for comparing the performance of different PBE solvers.

High‐order finite difference schemes for incompressible flows

AbstractThis paper presents a new high‐order approach to the numerical solution of the incompressible Stokes and Navier–Stokes equations. The class of schemes developed is based upon a velocity–pressure–pressure gradient formulation, which allows: (i) high‐order finite difference stencils to be applied on non‐staggered grids; (ii) high‐order pressure gradient approximations to be made using standard Padé schemes, and (iii) a variety of boundary conditions to be incorporated in a natural manner. Results are presented in detail for a selection of two‐dimensional steady‐state test problems, using the fourth‐order scheme to demonstrate the accuracy and the robustness of the proposed methods. Furthermore, extensions to higher orders and time‐dependent problems are illustrated, whereas the extension to three‐dimensional problems is also discussed. Copyright © 2010 John Wiley & Sons, Ltd.

Scalar Wave Equation Modeling with Time-Space Domain Dispersion-Relation-Based Staggered-Grid Finite-Difference Schemes

The staggered-grid finite-difference (SFD) method is widely used in numerical modeling of wave equations. Conventional SFD stencils for spatial derivatives are usually designed in the space domain. However, when they are used to solve wave equations, it becomes difficult to satisfy the dispersion relations exactly. Liu and Sen (2009c) proposed a new SFD scheme for one-dimensional (1D) scalar wave equation based on the time–space domain dispersion relation and plane wave theory, which is made to satisfy the exact dispersion relation. This new SFD scheme has greater accuracy and better stability than a conventional scheme under the same discretizations. In this paper, we develop this new SFD scheme further for numerical solution of 2D and 3D scalar wave equations. We demonstrate that the modeling accuracy is second order when the conventional 2 M -th-order space-domain SFD and the second order time-domain finite-difference stencils are directly used to solve the scalar wave equation. However, under the same discretization, our 1D scheme can reach 2 M -th-order accuracy and is always stable; 2D and 3D schemes can reach 2 M -th-order accuracy along 8 and 48 directions, respectively, and have better stability. The advantages of the new schemes are also demonstrated with dispersion analysis, stability analysis, and numerical modeling.

An estimation of the sensitivity of numerical error norm using adjoint model

AbstractA posteriori estimation of the numerical error sensitivity to the local truncation error is addressed using adjoint model endowed with the information on the error field. The numerical error is estimated from the solution of the linear tangent model (LTM) or from a Richardson extrapolation. The local truncation error used in the LTM is obtained by the action of a high‐order finite‐difference stencil on the field computed by the main (low‐order accuracy) algorithm. Copyright © 2009 John Wiley & Sons, Ltd.

Finite-difference stable stencils based on least-square quadric fitting

In this paper a procedure to obtain new finite difference stencils is given. We investigate the particular branch of cases where the order of the finite difference stencil is lower than the amount of grid points. New finite difference stencils are obtained and compared with traditional ones from literature. General properties and advantages of the proposed stencils are investigated, in particular the improvements obtained in stability when solving ill-posed problems.

A Family of Sixth-Order Compact Finite-Difference Schemes for the Three-Dimensional Poisson Equation

We derive a family of sixth-order compact finite-difference schemes for the three-dimensional Poisson's equation. As opposed to other research regarding higher-order compact difference schemes, our approach includes consideration of the discretization of the source function on a compact finite-difference stencil. The schemes derived approximate the solution to Poisson's equation on a compact stencil, and thus the schemes can be easily implemented and resulting linear systems are solved in a high-performance computing environment. The resulting discretization is a one-parameter family of finite-difference schemes which may be further optimized for accuracy and stability. Computational experiments are implemented which illustrate the theoretically demonstrated truncation errors.