Conventional Finite Difference Method Research Articles

An overset-grid finite-difference algorithm for seismic wavefield propagations modelling in the polar coordinate system with a complex free-surface topography

Summary Nowadays, finite-difference method has been widely applied to simulate seismic wavefield propagation in a local seismic model with a complex topography by utilizing curvilinear grids and traction image method in Cartesian coordinates system. For a global seismic model in polar coordinate system, it is still a challenge for the conventional finite-difference method to deal with the grid singularity at the center and the topographic free surface at the top. In this study, we develop a finite-difference method in polar coordinates for seismic wavefield propagation in the 2D global model. In the proposed finite-difference method, the overset-grid algorithm is used to handle the grid singularity at the center, and the curvilinear grid technique as well as the traction image method are applied to implement free-surface boundary condition on the complex topography. The proposed finite-difference method is validated in flat and topographic free-surface models by comparing synthetic waveforms with reference solutions. The seismic wavefield propagation in a realistic Mars profile is computed by the proposed finite-difference method. The proposed finite-difference method is an efficient and accurate method for seismic wave modeling in the polar coordinate system with a complex free-surface topography.

A coupled FD-SPH method for shock-structure interaction and dynamic fracture propagation modeling

Shock wave propagation and their damage to solid structures, which involve complex multiphase and multiphysics phenomena, are difficult problems to address. In this paper, a three-dimensional graphics processing unit (GPU)-accelerated finite difference-smoothed particle hydrodynamics (FD-SPH) method was developed for the prediction of strongly compressible fluid flow and strong fluid-structure interactions. The conventional mesh-based finite difference method was used to simulate shock wave propagation, while the smoothed particle hydrodynamics method was employed to predict the dynamic behavior of solid materials. In addition, the immersed boundary method (IBM) was implemented in the GPU-accelerated FD-SPH solver for the boundary treatment of the interface between solid and fluid materials. Four numerical cases, namely shock-cylinder obstacle interaction, elastic panel deformation induced by shock wave propagation, the response of stainless steel tubes under internal blast loading, and the damage of blast loading on reinforced concrete slab, were conducted for verification of the GPU-accelerated FD-SPH solver. The numerical results were compared against the available experimental data, which shows that the GPU-accelerated FD-SPH solver is capable of capturing the shock wave propagation and its damage to structures very well.

Numerical Simulation of Wave Equation Using the Finite Difference Method with Variable Refined Grid in Complex Near-surface Condition

The horizontal layered refined grid is adopted in the conventional finite difference method (FDM) of variable refined grid for the numerical simulation of wave equation. But, the method cannot well adapt to the characteristics of undulated topography and the variation of velocity structure in near-surface region. To solve this problem, a numerical simulation of wave equation using FDM with undulated variable refined grid is proposed in this paper. The grids are generated according to the undulated topography and the spatial distribution of velocity in the near-surface region from low-velocity to the high-velocity layer. For ensuring the accuracy while taking into account the computational efficiency, the FDM with variable coefficient in above grids is used to discretize the acoustic wave equation. The numerical evaluation and example show that the same accuracy as the method using fine grid can be obtained by the new method. However, its computational time is only 43.6% of the one spent by the fine grid. Furthermore, the new method can also stably adapt to the real complex loess plateau near-surface model in a working area of Ordos Basins.

Finite-difference method for modeling the surface wave propagation with surface topography in anisotropic-viscoelastic media

The growth of urbanization makes it increasingly important to accurately investigate underground spaces. Surface-wave-based inversion methods are becoming popular due to their cost-effectiveness and efficiency in reconstructing underground space. However, these methods may produce inaccurate results in the presence of media anisotropy or/and viscoelasticity. To investigate the influence of media anisotropy or viscoelasticity on surface waves, it is crucial to develop an efficient and accurate forward modeling method in anisotropic-viscoelastic (AV) media. The finite-difference (FD) method is a widely used forward modeling method in surface-wave-based inversion methods. However, implementing the free-surface boundary condition and representing non-flat topography can be challenging. Considering the challenges described above, this study proposes a simple and efficient FD method for modeling surface-wave propagation. The proposed method employs a parameter-modified strategy to implement the stress-free condition by modifying the model parameters near the (non-) flat free-surface boundary. A staircase discretization strategy is then used to represent the non-flat free surface. To suppress the staircase diffractions caused by the staircase discretization strategy, an independent wavefield superposition is adopted with modeling results of different parameter configurations. Compared with conventional FD methods, the proposed method can reduce computational costs while accurately simulating surface waves. Numerical experiments have demonstrated the accuracy and feasibility of the proposed method. The proposed method provides a theoretical basis for surface-wave-based inversion methods and contributes to the development of accurate underground space investigation.

Electro-Thermal Modelling by Novel Variational Methods: Racetrack Coil in Short-Circuit

The design of superconducting applications containing windings of superconducting wires or tapes requires electro-thermal quench modelling. In this article, we present a reliable numerical method based on a variational principle and we benchmark it to a conventional finite difference method that we implemented in C++. As benchmark problem, we consider a racetrack coil made of REBCO superconducting tape under short circuit, approximated as a DC voltage that appears at the initial time. Results show that both models agree with each other and analytical limits. Since both models take screening currents into account, they are promising for the design of magnets (especially fast-ramp magnets) and power applications, such as the stator windings of superconducting motors or generators.

The Method of Optical Paths for the Numerical Solution of Integrated Photonics Problems

A number of topical problems of integrated photonics are reduced to oblique incidence of radiation on a plane-parallel scatterer. For such problems, a method for integrating Maxwell’s equations along the direction of beam propagation is proposed. As a result, the original two-dimensional problem is reduced to a one-dimensional problem, and it is solved using recently proposed one-dimensional bicompact schemes. This significantly reduces the computational cost compared with the conventional two-dimensional finite difference and finite element methods. The proposed method is verified by solving test problems for which exact solutions are known.

Automatic differentiation for orbital-free density functional theory.

Differentiable programming has facilitated numerous methodological advances in scientific computing. Physics engines supporting automatic differentiation have simpler code, accelerating the development process and reducing the maintenance burden. Furthermore, fully differentiable simulation tools enable direct evaluation of challenging derivatives-including those directly related to properties measurable by experiment-that are conventionally computed with finite difference methods. Here, we investigate automatic differentiation in the context of orbital-free density functional theory (OFDFT) simulations of materials, introducing PROFESS-AD. Its automatic evaluation of properties derived from first derivatives, including functional potentials, forces, and stresses, facilitates the development and testing of new density functionals, while its direct evaluation of properties requiring higher-order derivatives, such as bulk moduli, elastic constants, and force constants, offers more concise implementations than conventional finite difference methods. For these reasons, PROFESS-AD serves as an excellent prototyping tool and provides new opportunities for OFDFT.

Efficient temporal high-order staggered-grid scheme with a dispersion-relation-preserving method for the scalar wave modeling

Staggered-grid finite-difference (FD) method is widely used to solve the wave equation for the numerical seismic modeling, and it is a key step of the advanced seismic imaging and inversion problem. However, the conventional FD method is prone to instability and dispersion error due to the insufficient approximation accuracy. In this work, we propose an efficient temporal high-order finite-difference (FD) scheme with the cross-rhombus stencil. Compared with the standard cross-rhombus method, the new method has less computational cost due to we simplify the FD scheme. Moreover, the dispersion relation of the new method is easy to obtain the dispersion-relation-preserving (DRP) FD coefficients, which can significantly alleviate the spatial and temporal dispersion errors. Dispersion and stability analyses indicate that the new scheme has better performance in seismic modeling than the conventional method, and numerical experiments also indicate that the new scheme can effectively mitigate dispersion error and improve the numerical accuracy.

Improved Finite Difference Technique via Adomian Polynomial to Solve the Coupled Drinfeld’s–Sokolov–Wilson System

This study presents a new algorithm for effectively solving the nonlinear coupled Drinfeld’s–Sokolov–Wilson (DSW) system using a hybrid explicit finite difference technique with the Adomian polynomial (EFD‐AP). The suggested approach addresses the problem of accurately solving the DSW system. Numerical results are obtained by comparing the exact solution with absolute and mean square errors using a test problem to assess the EFD‐AP accuracy against the exact solution and the conventional explicit finite difference (EFD) method. The results exhibit excellent agreement between the approximate and exact solutions at different time values, and the results showed that the proposed EFD‐AP method achieves superior accuracy and efficiency compared to the EFD method, which makes it a promising method for solving nonlinear partial differential systems of higher order.

Frequency-domain acoustic full waveform inversion with an embedded boundary method for irregular topography

In the implementation of full waveform inversion (FWI) to identify subsurface velocity distributions with land seismic data, which are often acquired in regions with irregular topography, wave equation-based modelling requires caution. In particular, when using the finite difference method (FDM), unwanted scattered waves are generated because irregular surfaces crossing a rectangular grid are discretized via a staircase approximation; hence, if the problems caused by this staircase approximation are disregarded, FDM-based FWI may fail due to the presence of undesirable wavefields. To resolve this problem, this study develops a 2D frequency-domain acoustic FWI technique using a 9-point FDM-based modelling scheme that includes an embedded boundary method (EBM). This study suggests a workflow for the whole EBM-based FWI process from the calculation of coefficients for the EBM-based 9-point FDM modelling to applying it to FWI for proper velocity updates. In numerical examples, using velocity models with a tilted surface and an arbitrarily fluctuating surface, we synthesize seismic data and verify the accuracy of EBM-based 9-point FDM modelling and its superiority over the conventional FDM by comparing it with wavefields derived from the spectral element method. Then, we show that our EBM-based FWI is able to estimate subsurface velocity distributions even though the model has irregular topography, which spoils the result of the conventional FWI.

An efficient symplectic stereo-modeling method for seismic inversion by using deep learning technique

AbstractThis paper proposes an efficient symplectic stereo-modeling (SSTEM) method for full waveform inversion (FWI) by using a deep learning technique. To solve the 2D acoustic equation, the SSTEM method uses a third-order optimal symplectic partitioned Runge–Kutta approach as a time-stepping method. An eighth-order stereo-modeling operator is used for spatial discretization. The SSTEM method is then expressed with a recurrent neural network (RNN). This is realized mainly because the time advancing format of the SSTEM method is similar to that of RNN, and they both use the information from the previous time step to obtain information from the current time step. With SSTEM as the forward modeling method, FWI is implemented using Tensorflow. The well-known adaptive moment estimation (Adam) optimizer and Nesterov adaptive moment estimation (Nadam) optimizer with mini-batch are used. The applicability of the developed code is also verified on GPUs. The numerical results show that the SSTEM method is more efficient and produces less numerical dispersion than the conventional finite-difference (FD) method when the same sampling rate in a wavelength is used. We compare several loss functions. The mean square (MSE) error and absolute (ABS) error loss functions are first tested. Another loss function that adds a physical differential operator to the original loss function is then considered. The FWI results show that this loss function has some improvements. Finally, we implement FWI on the complex Marmousi and SEG/EAGE models, and the inversion results demonstrate that the proposed method is suitable for seismic imaging in complex media.

An optimized finite-difference method to minimize numerical dispersion of acoustic wave propagation using a genetic algorithm

The finite-difference method (FDM) is one of the most popular methods for numerical simulation of wave propagation. The major challenge that we face in solving a partial differential equation for the wave propagation is numerical dispersion. We use an automated and optimized FDM using a genetic algorithm (GA) to optimally compute second-order spatial derivatives. In our method, the explicit finite-difference coefficients are calculated using the GA to minimize the dispersion (phase velocity) for all wavenumbers without using any specific window function. The amplitudes of the pseudospectral window (finite difference [FD] coefficients) are optimized by making the phase velocity close to the analytical solution at each wavenumber and stability close to that of the conventional FDM. Although FD coefficients in this method depend on velocity, grid spacing, and time step, less dispersive solutions can be achieved by computing suitable FD coefficients for varying cases. It takes only a few seconds of additional computation time depending on the order of approximation (e.g., approximately 20 s for the 12th order of approximation in a machine with 120 GB RAM and 2.0 GHz processor). Comparison of the numerical results from our algorithm with those from an existing pseudospectral method, the conventional FDM, and an existing time-space optimized FDM indicates that the proposed method is less dispersive for any order of approximation.

A New Stable, Explicit, Third‐Order Method for Diffusion‐Type Problems

AbstractThis paper reports on a novel explicit numerical method for the spatially discretized diffusion or heat equation. After discretizing the space variables as in conventional finite difference methods, this method does not use a finite difference approximation for the time derivatives, it instead combines constant‐neighbor and linear‐neighbor approximations, which decouple the ordinary differential equations, thus they can be solved analytically. In the obtained three‐stage method, the time step size appears in exponential form with negative coefficients in the final expression. This property guarantees unconditional stability, as it is shown using von Neumann stability analysis. It is also proved that the convergence of the method is third order in the time step size. After verification, by solving Fisher's and Huxley's equations, it is demonstrated that it works for nonlinear equations as well. The new algorithm is tested against widely used numerical solvers for cases where the media is strongly inhomogeneous. According to the results, the new method is significantly more effective than the traditional explicit or implicit methods, especially for extremely large stiff systems. It is believed that this new method is unique in the sense that it is the first unconditionally stable explicit method with third‐order convergence.

Multilayer heat equations: Application to finance

<p style='text-indent:20px;'>In this paper, we develop a Multilayer (ML) method for solving one-factor parabolic equations. Our approach provides a powerful alternative to the well-known finite difference and Monte Carlo methods. We discuss various advantages of this approach, which judiciously combines semi-analytical and numerical techniques and provides a fast and accurate way of finding solutions to the corresponding equations. To introduce the core of the method, we consider multilayer heat equations, known in physics for a relatively long time but never used when solving financial problems. Thus, we expand the analytic machinery of quantitative finance by augmenting it with the ML method. We demonstrate how one can solve various problems of mathematical finance by using our approach. Specifically, we develop efficient algorithms for pricing barrier options for time-dependent one-factor short-rate models, such as Black-Karasinski and Verhulst. Besides, we show how to solve the well-known Dupire equation quickly and accurately. Numerical examples confirm that our approach is considerably more efficient for solving the corresponding partial differential equations than the conventional finite difference method by being much faster and more accurate than the known alternatives.</p>

Rapid Surrogate Modeling of Electromagnetic Data in Frequency Domain Using Neural Operator

The efficiency of solving geophysical inverse problem largely relies on the efficiency of solving the corresponding forward problem. As for electromagnetic (EM) data forward modeling in frequency domain, the conventional numerical methods, e.g. finite difference method (FDM), discretize the governing equations resulting a large linear system which is usually expensive to solve. Meanwhile, for inversion iteration we normally do not need to solve the forward problem in high precision. Thus a rapid surrogate modeling approach which uses the neural network is promising for replacing the forward modeling module in the inversion scheme. Here we proposed an algorithm which uses the neural operator to solve the EM data modeling problem in the frequency domain. To develop a surrogate model for EM data forward problem, we introduce an extended Fourier neural operator (EFNO) that enables the calculation at least 100 times faster than the conventional FDM solver while maintaining good precision. Moreover, by adding a sub-network the proposed neural operator has good generalization which has the capacity of predicting solution at any site locations and frequencies. Due to the discretization-invariance of Fourier neural operator, the neural operator trained on coarse grids can easily transfer to fine grids with only retraining part of parameters, resulting in a super-resolution prediction capability. We test our proposed method with 2-D and 3-D magnetotelluric (MT) data modeling problems, demonstrating that the EFNO has great potentials for severing as a general rapid surrogate forward solver in EM data inversion scheme.

Wave-packet behaviors of the defocusing nonlinear Schrödinger equation based on the modified physics-informed neural networks.

Investigated in this paper is the defocusing nonlinear Schrödinger (NLS) equation, which is used for describing the wave-packet dynamics in certain weakly nonlinear media. With the physics-informed neural networks (PINNs), we modify the corresponding loss function in the existing literature and obtain two types of dark solitons, type-I and type-II solitons. It is demonstrated that the modified loss function presents higher-precision wave-packet behaviors based on fewer initial and boundary data. Taking type-I solitons into consideration, we find that when only a small fraction of initial and boundary data are given, the prediction accuracy of the wave packets will be increased one or two orders of magnitude at least if the modification term of the loss function is introduced. Furthermore, for the inverse problem, the modified loss function provides a better estimate of the nonlinear coefficient of the NLS equation based on fewer observed data of the wave packets. For type-II solitons, we compare the required data and predicted results of the PINNs with those of the conventional time-splitting finite difference (TSFD) method and reveal that achieving the same precision of the wave-packet behavior, the PINNs with the modified loss functions require only one tenth of the amount of the initial and boundary data of the TSFD method. Besides, both unmodified and modified loss functions are exploited for predicting the behaviors of Gaussian wave packets, and it is observed that the predicted result of the modified loss function agrees with the high-precision solution of the time-splitting Fourier pseudospectral method, whereas the unmodified loss function fails.

Simulating and imaging wavefield of irregular topography based on spectral element method

For the high-precision seismic imaging of complex structures on irregular topographies, this study proposes a Chebyshev spectral-element reverse time migration (CSE-RTM) method coupled with the implicit Newmark time integral method. The approach is first applied to a numerical example to consider factors affecting numerical accuracy, such as the time step and the pattern of subdivision of the computational domain. The results show that the smaller the time step is, the higher the numerical accuracy is and that numerical accuracy improves when the grid nodes tend to be distributed uniformly in the computational domain. The spectral element method combines the advantages of the stronger boundary-related adaptability of the finite element method with the high accuracy and fast convergence of the spectral method. This method is applied to reverse time migration imaging, and tests of the model show that its accuracy of imaging is higher than that of the conventional finite-difference reverse time migration method in case of complex structures in areas with irregular topography. To improve the efficiency of parallel computing of the CSE-RTM, an iterative algorithm is proposed to solve the hierarchical equilibrium sparse linear equations (HELPs) for it based on the idea of the degrees of cohesion of freedom and local relaxation. This algorithm has the advantage that the efficiency of parallel computation does not decrease with an increasing number of processors while the number of iterations required for convergence is controlled. Owing to the fast convergence of the spectral method on a sparse grid coupled with the parallel computational scheme of HELPs, the efficiency of calculation of the CSE-RTM is significantly improved.

A Set of New Stable, Explicit, Second Order Schemes for the Non-Stationary Heat Conduction Equation

This paper introduces a set of new fully explicit numerical algorithms to solve the spatially discretized heat or diffusion equation. After discretizing the space and the time variables according to conventional finite difference methods, these new methods do not approximate the time derivatives by finite differences, but use a combined two-stage constant-neighbour approximation to decouple the ordinary differential equations and solve them analytically. In the final expression for the new values of the variable, the time step size appears not in polynomial or rational, but in exponential form with negative coefficients, which can guarantee stability. The two-stage scheme contains a free parameter p and we analytically prove that the convergence is second order in the time step size for all values of p and the algorithm is unconditionally stable if p is at least 0.5, not only for the linear heat equation, but for the nonlinear Fisher’s equation as well. We compare the performance of the new methods with analytical and numerical solutions. The results suggest that the new algorithms can be significantly faster than the widely used explicit or implicit methods, particularly in the case of extremely large stiff systems.

3D Seismic-Wave Modeling with a Topographic Fluid–Solid Interface at the Sea Bottom by the Curvilinear-Grid Finite-Difference Method

ABSTRACTThe curvilinear-grid finite-difference method (FDM), which uses curvilinear coordinates to discretize the nonplanar interface geometry, is extended to simulate acoustic and seismic-wave propagation across the fluid–solid interface at the sea bottom. The coupled acoustic velocity-pressure and elastic velocity-stress formulation that governs wave propagation in seawater and solid earth is expressed in curvilinear coordinates. The formulation is solved on a collocated grid by alternative applications of forward and backward MacCormack finite difference within a fourth-order Runge–Kutta temporal integral scheme. The shape of a fluid–solid interface is discretized by a curvilinear grid to enable a good fit with the topographic interface. This good fit can obtain a higher numerical accuracy than the staircase approximation in the conventional FDM. The challenge is to correctly implement the fluid–solid interface condition, which involves the continuity of tractions and the normal component of the particle velocity, and the discontinuity (slipping) of the tangent component of the particle velocity. The fluid–solid interface condition is derived for curvilinear coordinates and explicitly implemented by a domain-decomposition technique, which splits a grid point on the fluid–solid interface into one grid point for the fluid wavefield and another one for the solid wavefield. Although the conventional FDM that uses effective media parameters near the fluid–solid interface to implicitly approach the boundary condition conflicts with the fluid–solid interface condition. We verify the curvilinear-grid FDM by conducting numerical simulations on several different models and compare the proposed numerical solutions with independent solutions that are calculated by the Luco-Apsel-Chen generalized reflection/transmission method and spectral-element method. Besides, the effects of a nonplanar fluid–solid interface and fluid layer on wavefield propagation are also investigated in a realistic seafloor bottom model. The proposed algorithm is a promising tool for wavefield propagation in heterogeneous media with a nonplanar fluid–solid interface.

Numerical Investigation of Progressive Slope Failure Induced by Sublevel Caving Mining Using the Finite Difference Method and Adaptive Local Remeshing

Slope failure induced by sublevel caving mining is a progressive process, resulting in the large deformation and displacement of rock masses in the slope. Numerical methods are widely used to investigate the above phenomenon. However, conventional numerical methods have difficulties when simulating the process of progressive slope failure. For example, the discrete element method (DEM) for block systems is computationally expensive and possibly fails for large-scale and complex slope models, while the finite difference method (FDM) has a mesh distortion problem when simulating progressive slope failure. To address the above problems, this paper presents a finite difference modeling method using the adaptive local remeshing technique (LREM) to investigate the progressive slope failure induced by sublevel caving mining. In the proposed LREM, (1) the zone of the distorted mesh is adaptively identified, and the landslide body is removed; (2) the updated mesh is regenerated by the local remeshing, and the physical field variables of the original computational model are transferred to the regenerated computational model. The novelty of the proposed method is that (1) compared with the DEM for block systems, the proposed LREM is capable of modeling the progressive slope failure in large-scale rock slopes; (2) the proposed method is able to address the problem of mesh distortion in conventional FDM modeling; and (3) compared with the errors induced by the frequent updating of the mesh of the entire model, the adaptive local remeshing technique effectively reduces calculation errors. To evaluate the effectiveness of the proposed LREM, it is first used to investigate the failure of a simplified slope induced by sublevel caving mining. Moreover, the proposed LREM is applied in a real case, i.e., to investigate the progressive slope failure induced by sublevel caving mining in Yanqianshan Iron Mine.