High-order Finite Difference Research Articles

On the role of global conservation property for finite difference schemes

High-order finite difference schemes are widely used in the field of numerical solutions of partial differential equations. Although the importance of their local conservation property is well recognized, especially for solving hyperbolic conservation laws, the importance of their global conservation property (GCP) is seldom considered. To show that in some cases the GCP is indeed important in improving performance of finite difference schemes, we make a comparative study of two fifth-order schemes with boundary closures satisfying the GCP or not. For some toy models considered in this short note, we show that the scheme satisfying the GCP does perform better than its counterpart without the GCP.

High-Order Bound-Preserving Finite Difference Methods for Multispecies and Multireaction Detonations

High-Order Bound-Preserving Finite Difference Methods for Multispecies and Multireaction Detonations

An immersed boundary fluid–structure interaction method for thin, highly compliant shell structures

A parallel computational method for simulating fluid–structure interaction problems involving large, geometrically nonlinear deformations of thin shell structures is presented and validated. A compressible Navier-Stokes solver utilizing a higher-order finite difference immersed boundary method is coupled with a geometrically nonlinear computational structural dynamics solver employing the mixed interpolation of tensorial components formulation for thin triangular shell elements. A weak fluid–structure coupling strategy is used to advance the numerical solution in time. The thin shell structures are represented in the fluid domain by a geometry mesh with a finite thickness at or below the size of the local grid spacing in the fluid domain. The methodologies for load and displacement transfer between the disparate geometry and structural meshes are detailed considering a parallel computing environment. The coupled method is validated for canonical simulation-based test cases and experimental fluid–structure interaction problems considering large deformations of thin shell structures, including a shock impinging on a cantilever plate, a fixed cylinder with a flexible trailing filament in channel flow, a thin, clamped plate in wall-bounded flow, and a flag waving in viscous crossflow. The FSI method is then demonstrated on a compliant circular sheet with a clamped center exposed to crossflow and finally applied to the inflation of a spacecraft disk-gap-band parachute inflating in supersonic flow conditions resembling the upper Martian atmosphere, where comparison with experimental data is provided.

An order-adaptive compact approximation Taylor method for systems of conservation laws

We present a new family of high-order shock-capturing finite difference numerical methods for systems of conservation laws. These methods, called Adaptive Compact Approximation Taylor (ACAT) schemes, use centered (2p+1)-point stencils, where p may take values in {1,2,…,P} according to a new family of smoothness indicators in the stencils. The methods are based on a combination of a robust first order scheme and the Compact Approximate Taylor (CAT) methods of order 2p-order, p=1,2,…,P so that they are first order accurate near discontinuities and have order 2p in smooth regions, where (2p+1) is the size of the biggest stencil in which large gradients are not detected. CAT methods, introduced in [3], are an extension to nonlinear problems of the Lax-Wendroff methods in which the Cauchy-Kovalesky (CK) procedure is circumvented following the strategy introduced in [22] that allows one to compute time derivatives in a recursive way using high-order centered differentiation formulas combined with Taylor expansions in time. The expression of ACAT methods for 1D and 2D systems of conservative laws is given and the performance is checked in a number of test cases for several linear and nonlinear systems of conservation laws, including Euler equations for gas dynamics.

Energy budget in decaying compressible MHD turbulence

Abstract

Structure preserving — Field directional splitting difference methods for nonlinear Schrödinger systems

A computational framework of high order conservative finite difference methods to approximate the solution of a general system of N coupled nonlinear Schrödinger equations (N-CNLS) is proposed. Exact conservation of the discrete analogues of the mass and the system’s Hamiltonian is achieved by decomposing the original system into a sequence of smaller nonlinear problems, associated to each component of the complex field, and a modified Crank–Nicolson time marching scheme appropriately designed for systems. For a particular model problem, we formally prove that a method, based on the standard second order difference formula, converges with order τ+h2; and, using the theory of composition method, schemes of order τ2+h2 and τ4+h2 are derived. The methodology can be easily extended to other high order finite difference formulas and composition methods. Conservation and accuracy are numerically validated.

Applying an advanced temporal and spatial high-order finite-difference stencil to 3D seismic wave modeling

Present finite-difference (FD) algorithms for modeling seismic wave propagation in 3D acoustic media are mainly based on the traditional orthogonal stencil and can achieve spatial high-order but temporal second-order accuracy. Therefore, these approaches may result in visible temporal dispersion and even instability for relatively large time sampling intervals. In this paper, we develop an advanced FD stencil with temporal and spatial arbitrary even-order accuracy to solve the 3D acoustic wave equation. The temporal derivative is approximated by an octahedron stencil with the operator length N together with the traditional second-order FD. Meanwhile, the spatial derivatives are calculated by the orthogonal stencil with the operator length M. The plane-wave theory is introduced to derive the time-space-domain dispersion relation and estimate the corresponding high-order FD coefficients by Taylor expansion. To achieve higher accuracy, we further optimize FD coefficients by applying least-squares (LS) to minimize the relative error of the time-space-domain dispersion relation. Dispersion and stability analyses show that our new schemes have higher modeling accuracy and propagation stability than the conventional method. Numerical examples demonstrate that the proposed FDs can generate greater temporal accuracy on the prerequisite of preserving satisfactory spatial accuracy. Our new temporal high-order FD stencil by incorporating larger time step, shorter operator lengths, LS-based coefficients and variable-length operators can be regarded as a kind of accurate and efficient tool for seismic wave modeling.

On the pseudospectral method and spectral accuracy

There are numerical accuracy problems related to the implementation of sharp internal interfaces in pseudospectral and finite-difference schemes. It is common practice to classify numerical errors due to the implementation of interfaces as being to some order in a Taylor expansion. An alternative approach is to classify these errors as being to some order in a Fourier expansion. The pseudospectral method does not provide spectral accuracy in inhomogeneous media. The numerical errors for the upper half of the frequency/wavenumber spectra of the propagating fields are not related to the implementation of the derivative operators but to aliasing effects coming from the multiplication of static material-parameter fields with the dynamic propagating fields. The pseudospectral method can only provide half-spectral accuracy. The same type of spatial aliasing errors is also present for finite-difference schemes. High-order finite differences can provide the same accuracy as the pseudospectral method if the staggered finite-difference derivative operators have a negligible error at four grid points per shortest wavelength and above. Smoothing of the material-parameter field leads to additional reduction in the error-free bandwidth of the propagating fields. Assuming that there is a maximum wavenumber up to which the spectrum of the smoothed model coincides with the implementation using a properly band-limited Heaviside step function, then there exists a local critical wavenumber for the propagating field equal to one half of the maximum wavenumber for the smoothed model. Harmonic averaging of material-parameter fields also results in wavenumber spectra where there is a maximum wavenumber above which the wavenumber spectrum deviates from an implementation with a band-limited Heaviside step function. The same one-half rule is also applicable in this case.

Uniformly convergent higher-order finite difference scheme for singularly perturbed parabolic problems with non-smooth data

Journal of Applied Mathematics and Computational Mechanics, Prace Naukowe Instytutu Matematyki i Informatyki, Politechnika Częstochowska, Scientific Research of the Institute of Mathematics and Computer Science, Czestochowa University of Technology

The design of new high-order group iterative method in the solution of two-dimensional fractional cable equation

In this article, a new high-order explicit group iterative scheme is developed for the solution of the two-dimensional time fractional cable equation. The proposed scheme is derived from the Crank-Nicolson (C-N) high-order compact finite difference method, where the Caputo discretization and C-N high-order approximations are used for the time fractional and space derivative respectively. The stability and convergence of the proposed scheme are also established. Finally, two numerical examples are presented to show the accuracy and efficiency of the scheme.

Very high-order Cartesian-grid finite difference method on arbitrary geometries

An arbitrary order finite difference method for curved boundary domains with Cartesian grid is proposed. The technique handles in a universal manner Dirichlet, Neumann or Robin conditions. We introduce the Reconstruction Off-site Data (ROD) method, that transfers in polynomial functions the information located on the physical boundary to the computational domain. Three major advantages are: (1) a simple description of the physical boundary with Robin condition using a collection of points; (2) no analytical expression (implicit or explicit) is required, particularly the ghost cell centroids projection are not needed; (3) we split up into two independent machineries the boundary treatment and the resolution of the interior problem, coupled by the ghost cell values. Numerical evidences based on the simple 2D convection-diffusion operators are presented to prove the capability of the method to reach at least the 6th-order with arbitrary smooth domains.

Variable-order optimal implicit finite-difference schemes for explicit time-marching solutions to wave equations

The time-space-domain finite-difference method (FDM) is widely used in forward modeling of wave equations. Conventional explicit FDMs (EFDMs) (with high order in space and the second order in time) would result in apparent temporal and spatial dispersion for high frequencies and large time steps. Moreover, the saturation effect of Taylor expansion seriously restricts the improvement of bandwidth coverage and efficiency of EFDMs. We have developed a variable-order optimization scheme of the implicit FDM (IFDM) to improve the computational efficiency and numerical accuracy of forward modeling. Then, we applied time-dispersion transforms (TDTs) to filter out the temporal dispersion generated by the second-order temporal approximations. Our method greatly alleviates the saturation effect of the high-order spatial finite-difference (FD) operators. Dispersion analysis indicates that the optimized coefficients of the IFDM based on the Remez algorithm can achieve the widest bandwidth (close to the Nyquist wavenumber), which corresponds to the shortest length of the spatial FD operator under a given error threshold. Our method has great potential to approach the highest spectral accuracy but with minimal increase in computational cost. Numerical experiments indicate that the combination of the variable-order optimization of IFDM and TDTs can significantly reduce numerical dispersion. Compared with traditional methods, our scheme is more advantageous for numerical simulation of large-scale geologic models because it has the least amount of calculation burden under the same accuracy requirements.

A Block–Interface Approach for High–Order Finite–Difference Simulations of Compressible Flows

The application of the high-order accurate schemes with multi-block domains is essential in problems with complex geometries. Primarily, accurate block-interface treatment is found to be of significant importance for precisely capturing discontinuities in such complex configurations. In the current study, a conservative and accurate multi-block strategy is proposed and implemented for a high-order compact finite-difference solver. For numerical discretization, the Beam-Warming linearization scheme is used and further extended for three-dimensional problems. Moreover, the fourth-order compact finite-difference scheme is employed for spatial discretization. The capability of the high-order multi-block approach is then evaluated for the onedimensional flow inside a Shubin nozzle, two-dimensional flow over a circular bump, and three-dimensional flow around a NACA 0012 airfoil. The results showed a reasonable agreement with the available exact solutions and simulation results in the literature. Further, the proposed block-interface treatment performed quite well in capturing shock waves, even in situations that the location of the shock coincides with block interfaces.

Large-Eddy Simulation of Compressible Flow in Asymmetric and Symmetric Planar Nozzles

Large-eddy simulations of supersonic, turbulent flow in asymmetric and symmetric planar nozzles are carried out using high-order finite difference schemes and a large-eddy simulation approach based on explicit filtering. The inflow to the nozzles is from supersonic, fully developed turbulent channel flows at and 360 with centerline Mach number 1.2. The combined effects of expansion and acceleration of the flow in this geometry result in reductions of mean pressure, density, and temperature as well as a reduction of Reynolds stresses. The asymmetry of the nozzle leads to notable differences in the aforementioned effects near the straight and the curved walls. It is shown that the decay in Reynolds stresses is more near the curved wall than near the straight wall of asymmetric nozzle, and it is nearly the same near both walls of the symmetric nozzle. Effects of longitudinal streamline curvature are found to be localized near the curved wall of the asymmetric nozzle and are quantified. The aforementioned effects of acceleration, expansion, and streamline curvature are found to be similar in the flow cases with two different Reynolds numbers, although with increasing Reynolds number, subtle differences in these effects are found.

Properties of High-Order Finite Difference Schemes and Idealized Numerical Testing

Construction of high-order difference schemes based on Taylor series expansion has long been a hot topic in computational mathematics, while its application in comprehensive weather models is still very rare. Here, the properties of high-order finite difference schemes are studied based on idealized numerical testing, for the purpose of their application in the Global/Regional Assimilation and Prediction System (GRAPES) model. It is found that the pros and cons due to grid staggering choices diminish with higher-order schemes based on linearized analysis of the one-dimensional gravity wave equation. The improvement of higher-order difference schemes is still obvious for the mesh with smooth varied grid distance. The results of discontinuous square wave testing also exhibits the superiority of high-order schemes. For a model grid with severe non-uniformity and non-orthogonality, the advantage of high-order difference schemes is inapparent, as shown by the results of two-dimensional idealized advection tests under a terrain-following coordinate. In addition, the increase in computational expense caused by high-order schemes can be avoided by the precondition technique used in the GRAPES model. In general, a high-order finite difference scheme is a preferable choice for the tropical regional GRAPES model with a quasi-uniform and quasi-orthogonal grid mesh.

Diffuse interface approach to modeling wavefields in a saturated porous medium

We present a model and numerical technique method for simulation a small-amplitude wave propagation in the fluid-saturated porous medium in regions with different porosity, including domains with pure fluid and pure solid. The governing equations are derived from the general hyperbolic thermodynamically compatible model of compressible fluid flow in a deformed porous medium. The resulting linear partial differential equations (PDE) form the symmetric hyperbolic system written for velocities, relative velocities, pressure and shear stress, which ensures application of a high order finite difference scheme on a staggered grid. A diffuse interface approach is applied to simulate wavefields in regions with interfaces between media with different porosity, including pure phases. The latter allows computations with the use of a single PDE system on rectangular grids without cumbersome interface tracking. Numerical experiments confirming the efficiency of the proposed approach are presented.

A high-order compact finite difference method on nonuniform time meshes for variable coefficient reaction–subdiffusion problems with a weak initial singularity

A high-order compact finite difference method on nonuniform time meshes is proposed for solving a class of variable coefficient reaction–subdiffusion problems. The solution of such a problem in general has a typical weak singularity at the initial time. Alikhanov’s high-order approximation on a uniform time mesh for the Caputo time fractional derivative is generalised to a class of nonuniform time meshes, and a fourth-order compact finite difference scheme is used for approximating the spatial variable coefficient differential operator. A full theoretical analysis of the stability and convergence of the method is given for the general case of the variable coefficients by developing an analysis technique different from the one for the constant coefficient problem. Taking the weak initial singularity of the solution into account, a sharp error estimate in the discrete $$L^{2}$$ -norm is obtained. It is shown that the proposed method attains the temporal optimal second-order convergence provided a proper mesh parameter is employed. Numerical results demonstrate the sharpness of the theoretical error analysis result.

A class of compact finite difference schemes for solving the 2D and 3D Burgers’ equations

In this paper, a class of two-level high-order compact finite difference implicit schemes are proposed for solving the Burgers’ equations. Firstly, based on the fourth-order compact finite difference schemes for spatial derivatives and the truncation error remainder correction method for temporal derivative, the high-order compact (HOC) difference method is introduced for solving the one-dimensional (1D) Burgers’ equation. At the same time, the stability of the scheme is analyzed by using the Fourier analysis method. Because only three grid points are involved in each time level. The Thomas algorithm can be directly used to solve the tridiagonal linear system. Then, this method is extended to solve the two-dimensional (2D) and three-dimensional (3D) coupled Burgers’ equations. Finally, numerical experiments are conducted to verify the accuracy and the reliability of the present schemes.

High-order accurate entropy stable finite difference schemes for the shallow water magnetohydrodynamics

This paper develops the high-order accurate entropy stable (ES) finite difference schemes for the shallow water magnetohydrodynamic (SWMHD) equations. They are built on the numerical approximation of the modified SWMHD equations with the Janhunen source term. First, the second-order accurate well-balanced semi-discrete entropy conservative (EC) schemes are constructed, satisfying the entropy identity for the given convex entropy function and preserving the steady states of the lake at rest (with zero magnetic field). The key is to match both discretizations for the fluxes and the non-flat river bed bottom and Janhunen source terms, and to find the affordable EC fluxes of the second-order EC schemes. Next, by using the second-order EC schemes as building block, high-order accurate well-balanced semi-discrete EC schemes are proposed. Then, the high-order accurate well-balanced semi-discrete ES schemes are derived by adding a suitable dissipation term to the EC scheme with the WENO reconstruction of the scaled entropy variables in order to suppress the numerical oscillations of the EC schemes. After that, the semi-discrete schemes are integrated in time by using the high-order strong stability preserving explicit Runge-Kutta schemes to obtain the fully-discrete high-order well-balanced schemes. The ES property of the Lax-Friedrichs flux is also proved and then the positivity-preserving ES schemes are studied by using the positivity-preserving flux limiters. Finally, extensive numerical tests are conducted to validate the accuracy, the well-balanced, ES and positivity-preserving properties, and the ability to capture discontinuities of our schemes.

High order nonlinear filter methods for subsonic turbulence simulation with stochastic forcing

Numerical simulations of forced turbulence in compressible fluids are challenging due to the multi-scale nature of the problem and conflicting requirements for numerical methods to accurately resolve the small scales and, at the same time, to handle shock waves and other discontinuities without generating spurious oscillations. Minimizing nonlinear instability and aliasing error while maintaining high accuracy before the simulation reaches the statistically stationary stage is computationally intensive. The goal of this work is to employ an efficient class of high order finite difference nonlinear filter methods for subsonic turbulence simulation with stochastic forcing. The 3D Euler equations for subsonic turbulence with temporally varying stochastic forcing at rms Mach numbers up to 0.6 are numerically solved using the Strang operator splitting of the homogeneous part of the Euler equations and the forcing terms. It was shown in Yee et al. (2013) that the Strang operator splitting is more stable than solving the full Euler equation with the source term included. The spatially seventh-order nonlinear filter methods with adaptive dissipation control developed by Yee & Sjögreen (2007, 2011) are used to solve the homogeneous system and an ODE solver is used to solve the forcing source terms. The nonlinear filter method includes a full time step of a spatially eighth-order central base method, using a third-order TVD Runge-Kutta time integration. The solution computed with the central base method is then nonlinearly filtered by an adaptive flow sensor and the dissipative portion of a seventh-order WENO with the Roe Riemann solver. In order to improve nonlinear stability of the base method without added numerical dissipation, the central base method discretizes the skew-symmetric split form of the inviscid flux derivatives. Both Ducros et al. and Kennedy-Gruber skew-symmetric split forms are tested. Numerical stability, computational efficiency, and effective spectral bandwidth of the nonlinear filter schemes are compared with those of second-order TVD and fifth- and seventh-order WENO methods. It is shown that the nonlinear filter method for this application is substantially more efficient, accurate and yields a superior spectral bandwidth compared to the standard TVD and WENO methods. The nonlinear filter method also demonstrates robust long-time integration for moderately compressible, statistically stationary turbulence with large-scale solenoidal forcing, including small-scale quantities such as enstrophy and mean-square dilatation.