High-order Finite-difference Operators Research Articles

Weak Versus Strong Wall Boundary Conditions for the Incompressible Navier-Stokes Equations

The pressure-velocity formulation of the incompressible Navier-Stokes equations is solved using high-order finite difference operators satisfying a summation-by-parts property. Two methods for imposing Dirichlet boundary conditions (one strong and one weak) are presented and proven stable using the energy method. Additionally, novel diagonal-norm second-derivative finite difference operators are derived with highly improved boundary accuracy. Accuracy and convergence measurements are presented and verified against theoretical expectations. Numerical experiments also show that subtle effects close to solid walls are more efficiently captured with strong boundary condition imposition methods rather than weak (less degrees of freedom required).

CLAIRE: Constrained Large Deformation Diffeomorphic Image Registration on Parallel Computing Architectures.

SummaryCLAIRE (Mang & Biros, 2019) is a computational framework for Constrained LArge deformation diffeomorphic Image REgistration (Mang et al., 2019). It supports highly-optimized, parallel computational kernels for (multi-node) CPU (Gholami et al., 2017; Mang et al., 2019; Mang & Biros, 2016) and (multi-node multi-)GPU architectures (Brunn et al., 2020, 2021). CLAIRE uses MPI for distributed-memory parallelism and can be scaled up to thousands of cores (Mang et al., 2019; Mang & Biros, 2016) and GPU devices (Brunn et al., 2020). The multi-GPU implementation uses device direct communication. The computational kernels are interpolation for semi-Lagrangian time integration, and a mixture of high-order finite difference operators and Fast-Fourier-Transforms (FFTs) for differentiation. CLAIRE uses a Newton–Krylov solver for numerical optimization (Mang & Biros, 2015, 2017). It features various schemes for regularization of the control problem (Mang & Biros, 2016) and different similarity measures. CLAIRE implements different preconditioners for the reduced space Hessian (Brunn et al., 2020; Mang et al., 2019) to optimize computational throughput and enable fast convergence. It uses PETSc (Balay et al., n.d.) for scalable and efficient linear algebra operations and solvers and TAO (Balay et al., n.d.; Munson et al., 2015) for numerical optimization. CLAIRE can be downloaded at https://github.com/andreasmang/claire.

Spectral properties of the incompressible Navier-Stokes equations

The influence of different boundary conditions on the spectral properties of the incompressible Navier-Stokes equations is investigated. By using the Fourier-Laplace transform technique, we determine the spectra, extract the decay rate in time, and investigate the dispersion relation. In contrast to an infinite domain, where only diffusion affects the convergence, we show that also the propagation speed influence the rate of convergence to steady state for a bounded domain. Once the continuous equations are analyzed, we discretize using high-order finite-difference operators on summation-by-parts form and demonstrate that the continuous analysis carries over to the discrete setting. The theoretical results are verified by numerical experiments, where we highlight the necessity of high accuracy for a correct description of time-dependent phenomena.

Multi-Node Multi-GPU Diffeomorphic Image Registration for Large-Scale Imaging Problems.

We present a Gauss-Newton-Krylov solver for large deformation diffeomorphic image registration. We extend the publicly available CLAIRE library to multi-node multi-graphics processing unit (GPUs) systems and introduce novel algorithmic modifications that significantly improve performance. Our contributions comprise (i) a new preconditioner for the reduced-space Gauss-Newton Hessian system, (ii) a highly-optimized multi-node multi-GPU implementation exploiting device direct communication for the main computational kernels (interpolation, high-order finite difference operators and Fast-Fourier-Transform), and (iii) a comparison with state-of-the-art CPU and GPU implementations. We solve a 2563-resolution image registration problem in five seconds on a single NVIDIA Tesla V100, with a performance speedup of 70% compared to the state-of-the-art. In our largest run, we register 20483 resolution images (25 B unknowns; approximately 152× larger than the largest problem solved in state-of-the-art GPU implementations) on 64 nodes with 256 GPUs on TACC's Longhorn system.

Performance and stability of the double absorbing boundary method for acoustic-wave propagation

The double absorbing boundary (DAB) is a novel extension to the family of high-order absorbing boundary condition operators. It uses auxiliary variables in a boundary layer to set up cancellation waves that suppress wavefield energy at the computational-domain boundary. In contrast to the perfectly matched layer (PML), the DAB makes no assumptions about the incoming wavefield and can be implemented with a boundary layer as thin as three computational grid-point cells. Our implementation incorporates the DAB into the boundary cell layer of high-order finite-difference (FD) techniques, thus avoiding the need to specify a padding region within the computational domain. We tested the DAB by propagating acoustic waves through homogeneous and heterogeneous 3D earth models. Measurements of the spectral response of energy reflected from the DAB indicate that it reflects approximately 10–15 dB less energy for heterogeneous models than a convolutional PML of the same computational memory complexity. The same measurements also indicate that a DAB boundary layer implemented with second-order FD operators couples well with higher-order FD operators in the computational domain. Long-term stability tests find that the DAB and CPML methods are stable for the acoustic-wave equation. The DAB has promise as a robust and memory-efficient absorbing boundary for 3D seismic imaging and inversion applications as well as other wave-equation applications in applied physics.

Practical inlet boundary conditions for internal flow calculations

To impose boundary conditions, data at the boundaries must be known, and consequently measurements of the imposed quantities must be available. In this paper, we consider the two most commonly used inflow boundary conditions with available data for internal flow calculations: the specification of the total temperature and total pressure. We use the energy method to prove that the specification of the total temperature and the total pressure together with the tangential velocity at an inflow boundary lead to well-posedness for the linearized compressible Euler equations. Next, these equations are discretized in space using high-order finite-difference operators on summation-by-parts form, and the boundary conditions are weakly imposed. The resulting numerical scheme is proven to be stable and the implementation of the corresponding nonlinear scheme is verified with the method of manufactured solutions. We also derive the spectrum for the continuous and discrete problems and show how to predict the convergence rate to steady state. Finally, nonlinear steady-state computations are performed, and they confirm the predicted convergence rates.

A Flexible Boundary Procedure for Hyperbolic Problems: Multiple Penalty Terms Applied in a Domain

AbstractA new weak boundary procedure for hyperbolic problems is presented. We consider high order finite difference operators of summation-by-parts form with weak boundary conditions and generalize that technique. The new boundary procedure is applied near boundaries in an extended domain where data is known. We show how to raise the order of accuracy of the scheme, how to modify the spectrum of the resulting operator and how to construct non-reflecting properties at the boundaries. The new boundary procedure is cheap, easy to implement and suitable for all numerical methods, not only finite difference methods, that employ weak boundary conditions. Numerical results that corroborate the analysis are presented.

An intrusive hybrid method for discontinuous two-phase flow under uncertainty

An intrusive stochastic projection method for two-phase time-dependent flow subject to uncertainty is presented. Numerical experiments are carried out assuming uncertainty in the location of the physical interface separating the two phases, but the framework generalizes to uncertainty with known distribution in other input data. Uncertainty is represented through a truncated multiwavelet expansion.We assume that the discontinuous features of the solution are restricted to computational subdomains and use a high-order method for the smooth regions coupled weakly through interfaces with a robust shock capturing method for the non-smooth regions.The discretization of the non-smooth region is based on a generalization of the HLL flux, and have many properties in common with its deterministic counterpart. It is simple and robust, and captures the statistics of the shock. The discretization of the smooth region is carried out with high-order finite-difference operators satisfying a summation-by-parts property.

High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains

Nonlinear entropy stability is used to derive provably stable high-order finite difference operators including boundary closure stencils, for the compressible Navier–Stokes equations. A comparison technique is used to derive a new Entropy Stable Weighted Essentially Non-Oscillatory (SSWENO) finite difference method, appropriate for simulations of problems with shocks. Viscous terms are approximated using conservative, entropy stable, narrow-stencil finite difference operators. The efficacy of the new discrete operators is demonstrated using both smooth and discontinuous test cases.

Staggered-Grid Finite Difference Method with Variable-Order Accuracy for Porous Media

The numerical modeling of wave field in porous media generally requires more computation time than that of acoustic or elastic media. Usually used finite difference methods adopt finite difference operators with fixed-order accuracy to calculate space derivatives for a heterogeneous medium. A finite difference scheme with variable-order accuracy for acoustic wave equation has been proposed to reduce the computation time. In this paper, we develop this scheme for wave equations in porous media based on dispersion relation with high-order staggered-grid finite difference (SFD) method. High-order finite difference operators are adopted for low-velocity regions, and low-order finite difference operators are adopted for high-velocity regions. Dispersion analysis and modeling results demonstrate that the proposed SFD method can decrease computational costs without reducing accuracy.

An improved vacuum formulation for 2D finite-difference modeling of Rayleigh waves including surface topography and internal discontinuities

Rayleigh waves are generated along the free surface and their propagation can be strongly influenced by surface topography. Modeling of Rayleigh waves in the near surface in the presence of topography is fundamental to the study of surface waves in environmental and engineering geophysics. For simulation of Rayleigh waves, the traction-free boundary condition needs to be satisfied on the free surface. A vacuum formulation naturally incorporates surface topography in finite-difference (FD) modeling by treating the surface grid nodes as the internal grid nodes. However, the conventional vacuum formulation does not completely fulfill the free-surface boundary condition and becomes unstable for modeling using high-order FD operators. We developed a stable vacuum formulation that fully satisfies the free-surface boundary condition by choosing an appropriate combination of the staggered-grid form and a parameter-averaging scheme. The elastic parameters on the topographic free surface are updated with exactly the same treatment as internal grid nodes. The improved vacuum formulation can accurately and stably simulate Rayleigh waves along the topographic surface for homogeneous and heterogeneous elastic models with high Poisson’s ratios ([Formula: see text]). This method requires fewer grid points per wavelength than the stress-image-based methods. Internal discontinuities in a model can be handled without modification of the algorithm. Only minor changes are required to implement the improved vacuum formulation in existing 2D FD modeling codes.

Finite difference order doubling in two dimensions

An order doubling process previously used to obtain eighth-order eigenvalues from the fourth-order Numerov method is applied to the perturbed oscillator in two dimensions. A simple method of obtaining high order finite difference operators is reported and an odd parity boundary condition is found to be effective in facilitating the smooth operation of the order doubling process.

Towards High Order Numerical Simulation of Aeolian Tones

AbstractStrictly stable high order finite difference operators have been applied to the compressible Navier‐Stokes equations in perturbation form for low Mach number computational aeroacoustics. Aeolian tones generated by vortex shedding from a circular cylinder have been simulated. (© 2005 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)

Differential quadrature finite difference method for structural mechanics problems

AbstractThe differential quadrature finite difference method (DQFDM) has been proposed by the author. The finite difference operators are derived by the differential quadrature (DQ). They can be obtained by using the weighting coefficients for DQ discretizations. The derivation is straight and easy. By using different orders or the same order but different grid DQ discretizations for the same derivative or partial derivative, various finite difference operators for the same differential or partial differential operator can be obtained. Finite difference operators for unequally spaced and irregular grids can also be generated through the use of generic differential quadrature (GDQ). The derivation of higher‐order finite difference operators is also easy. By adopting the same order of approximation for all mathematical terms existing in the problem to be solved, excellent convergence can be obtained due to the consistent approximation.The DQFDM is effective for solving structural mechanics problems. The numerical simulations for solving anisotropic nonuniform plate problems and two‐dimensional plane elasticity problems are carried out. Numerical results are presented. They demonstrate the DQFDM. Copyright © 2001 John Wiley Sons, Ltd.

Solution of anisotropic nonuniform plate problems by the differential quadrature finite difference method

The differential quadrature finite difference method (DQFDM) has been proposed by the author. The finite difference operators are derived by the differential quadrature (DQ). They can be obtained by using the weighting coefficients for DQ discretizations. The derivation is straight and easy. By using different orders or the same order but different grid DQ discretizations for the same derivative or partial derivative, various finite difference operators for the same differential or partial differential operator can be obtained. Finite difference operators for unequally spaced and irregular grids can also be generated through the use of generic differential quadrature (GDQ). The derivation of higher order finite difference operators is also easy. The DQFDM is used to solve anisotropic nonuniform plate problems. Numerical results are presented. They demonstrate the DQFDM.

Finite difference discretizations by differential quadrature techniques

A differential quadrature finite difference method (DQFDM) is proposed. The finite difference operators are derived by the differential quadrature (DQ). They can be obtained by using the weighting coefficients for DQ discretizations. The derivation is straightforward. By using different orders or the same order but different grid DQ discretizations for the same derivative or partial derivative, various finite difference operators for the same differential or partial differential operator can be obtained. Finite difference operators for unequally spaced and irregular grids can also be generated through the use of generic differential quadrature (GDQ). The derivation of higher order finite difference operators is also easy.

A Stable and Conservative Interface Treatment of Arbitrary Spatial Accuracy

Stable and accurate interface conditions based on the SAT penalty method are derived for the linear advection–diffusion equation. The conditions are functionally independent of the spatial order of accuracy and rely only on the form of the discrete operator. We focus on high-order finite-difference operators that satisfy the summation-by-parts (SBP) property. We prove that stability is a natural consequence of the SBP operators used in conjunction with the new, penalty type, boundary conditions. In addition, we show that the interface treatments are conservative. The issue of the order of accuracy of the interface boundary conditions is clarified. New finite-difference operators of spatial accuracy up to sixth order are constructed which satisfy the SBP property. These finite-difference operators are shown to admit design accuracy (pth-order global accuracy) when (p−1)th-order stencil closures are used near the boundaries, if the physical boundary conditions and interface conditions are implemented to at leastpth-order accuracy. Stability and accuracy are demonstrated on the nonlinear Burgers' equation for a 12-subdomain problem with randomly distributed interfaces.

Transient fluid-structure interaction in naval applications using the retarded potential method

Submerged structures subjected to underwater shock remain a difficult and costly problem to analyze. In this work, the retarded potential (RP) method, which is a boundary element technique, is employed to discretize the fluid domain. The retarded potential program RPF is coupled to the ABAQUS non-linear finite element program to form ABAQUS/RPF. The Figueiredo method is used to convert the standard RP differential-delay equations for the fluid to simply delay equations, which are inherently more stable. ABAQUS/RPF is applied to both a rigid and elastic sphere. The Figueiredo approach is seen to extend significantly the stability range of the standard RP by a factor of approximately 10–20. Integration of the constant-pressure RP fluid element is shown to be very important to achieve this level of accuracy. Efforts to improve the Figueiredo approach using higher-order finite difference operators appear to have no effect on accuracy or stability. Stability of the coupled system is achieved by simply not extrapolating the structural accelerations.

Numerical studies of hyperbolic IBVP with high-order finite difference operators satisfying a summation by parts rule

Numerical studies on hyperbolic initial-boundary value problems (IBVP) have been performed using high-order difference operators satisfying a summation by parts rule. To assure that the numerical solution is strictly stable two recently developed methods by Carpenter et al. (1994) and Olsson (1995) to implement the analytic boundary conditions without destroying the summation by parts rule have been used. Theoretical and numerical results show that the numerical methods presented here are strictly stable and have a convergence rate that agrees well with the theory of Gustafsson (1975, 1981).

The application of a two-dimensional higher-order theory for the analysis of a thick elastic plate

The application of approximate equations of a two-dimensional theory for the analysis of stress and displacement distributions of a thick elastic plate with fixed boundaries is presented. By using the method of power series expansion of displacement components, a set of fundamental equations of a two-dimensional higher-order plate theory is derived through the principle of virtual work. In order to assure the accuracy of the present theory, the governing equations of several sets of truncated approximate theories are solved numerically by the finite difference method. Higher-order finite difference operators with error terms of o( h 4) are applied in order to validate the numerical accuracy of results. The convergence properties of the present numerical solutions have been shown to be accurate for the expanded displacement components with respect not only to the number of finite difference mesh intervals but also to the order of approximate theories. Stress components are determined by satisfying the equilibrium equations of a three-dimensional elastic continuum and the stress boundary conditions on the upper and lower surfaces of a plate. A direct comparison is also made with the numerical results of the integral method obtained by using the solutions of a three-dimensional elastic continuum.