High-order Stencil Research Articles

Line-Based High-Order Hybrid Scheme on Unstructured Grids for Compressible Flows

This study explores the use of line-based finite difference techniques on unstructured grids employing high-order stencils for practical aerodynamic applications. The primary focus is on finite-difference explicit and compact spatial discretization up to sixth order. Problems studied include compressible flows containing discontinuities, which necessitated the use of a hybrid finite difference–finite volume approach, where the former is shown to be effective and accurate for smooth regions, while the latter is robust in capturing discontinuities and shocks. To differentiate between discontinuous and smooth zones, two distinct shock indicators are used, one based on pressure gradients and the other on the divergence of velocity. To achieve faster convergence, spatial discretization is applied with line-based alternating direction implicit (ADI) time integration. A key focus of this work is to demonstrate the proposed methodology in flow scenarios relevant to aircraft applications that also test various aspects of the methodology. Consequently, the cases studied are isentropic vortex convection, standing normal shock, transonic NACA0012 airfoil, supersonic cylinder, unsteady flow past tandem circular cylinders, and supersonic flow past a forward-facing step. Comparisons were also done against an existing line-based finite-volume technique for unstructured grids that employed Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) and fifth-order Weighted Essentially NonOscillatory (WENO5) spatial reconstruction schemes. It is observed that the hybrid finite difference–finite volume approach performed well with excellent agreement against available results, sometimes with far fewer degrees of freedom. Furthermore, the implicit ADI scheme is effective over unstructured grids and provides a three-time speed-up over the traditional Runge–Kutta fourth-order scheme.

Rotation measure structure functions with higher-order stencils as a probe of small-scale magnetic fluctuations and its application to the Small and Large Magellanic Clouds

ABSTRACTMagnetic fields and turbulence are important components of the interstellar medium (ISM) of star-forming galaxies. It is challenging to measure the properties of the small-scale ISM magnetic fields (magnetic fields at scales smaller than the turbulence driving scale). Using numerical simulations, we demonstrate how the second-order rotation measure (RM, which depends on thermal electron density, ne, and magnetic field, b) structure function can probe the properties of small-scale b. We then apply our results to observations of the Small and Large Magellanic Clouds (SMC and LMC). First, using Gaussian random b, we show that the characteristic scale, where the RM structure function flattens is approximately equal to the correlation length of b. We also show that computing the RM structure function with a higher-order stencil (more than the commonly-used two-point stencil) is necessary to accurately estimate the slope of the structure function. Then, using Gaussian random b and lognormal ne with known power spectra, we derive an empirical relationship between the slope of the power spectrum of b, ne, and RM. We apply these results to the SMC and LMC and estimate the following properties of small-scale b: correlation length (160 ± 21 pc for the SMC and 87 ± 17 pc for the LMC), strength (14 ± 2 $\mu {\rm G}$ for the SMC and 15 ± 3 $\mu {\rm G}$ for the LMC), and slope of the magnetic power spectrum (−1.3 ± 0.4 for the SMC and −1.6 ± 0.1 for the LMC). We also find that ne is practically constant over the estimated b correlation scales.

Generating surface-offset common-image gathers with backward wavefield synthesis

Surface-offset common-image gathers (CIGs) are an important data format for seismic velocity analysis. However, the reverse time migration (RTM) method, which is wavefield propagation based, does not directly produce surface-offset CIGs because it propagates waves from all the offsets together. Here, we implement the generation of surface-offset CIGs by synthesizing the backward wavefields using the forward wavefields at locations where source and receiver locations overlap. This method produces surface-offset CIGs with high accuracy and low cost when compared with those generated by backpropagating each trace separately. We adopt nonnegative least-squares filters for sparse linear deconvolution. The computational effectiveness of the proposed method increases for higher dimensions, higher-order stencils, and more complicated wave equations. The proposed method works stably on a realistic towed-streamer acquisition system with moderate geometric positioning errors between the locations of airguns and hydrophones.

Scalable communication for high-order stencil computations using CUDA-aware MPI

Modern compute nodes in high-performance computing provide a tremendous level of parallelism and processing power. However, as arithmetic performance has been observed to increase at a faster rate relative to memory and network bandwidths, optimizing data movement has become critical for achieving strong scaling in many communication-heavy applications. This performance gap has been further accentuated with the introduction of graphics processing units, which can provide by multiple factors higher throughput in data-parallel tasks than central processing units. In this work, we explore the computational aspects of iterative stencil loops and implement a generic communication scheme using CUDA-aware MPI, which we use to accelerate magnetohydrodynamics simulations based on high-order finite differences and third-order Runge–Kutta integration. We put particular focus on improving intra-node locality of workloads. Our GPU implementation scales strongly from one to 64 devices at 50%–87% of the expected efficiency based on a theoretical performance model. Compared with a multi-core CPU solver, our implementation exhibits 20–60× speedup and 9–12× improved energy efficiency in compute-bound benchmarks on 16 nodes.

Finite differences for higher order derivatives of low resolution data

The viability of finite difference stencils for higher-order derivatives of low-resolution data, encountered in various applications, is examined. This study’s illustrative example involves the evaluation of the fourth derivative of the digitized free surface radius obtained from pixelated images. Procedures to obtain an optimal approximation to the derivative are made from estimates of truncation and roundoff error, with emphasis on the analysis of roundoff error. A method of successive approximations to evaluate the optimal grid spacing is presented. Higher-order stencils allow for larger optimal grid spacing, thereby reducing the roundoff error that would otherwise dominate. Hence, higher-order stencils are effective for low precision data.

Modeling 3D acoustic-wave propagation using modified cuboid-based staggered-grid finite-difference methods with temporal and spatial high-order accuracy

To improve the modeling accuracy and adaptability of traditional temporal second-order staggered-grid finite-difference (SFD) methods for 3D acoustic-wave modeling, we propose a modified time-space-domain temporal and spatial high-order SFD stencil on a cuboid grid. The grid nodes on a double-pyramid stencil and the standard orthogonality stencil are used to approximate temporal and spatial derivatives. This stencil can adopt different grid spacing in each spatial axis, and thus it is more flexible than the existing one with the same grid spacing. Based on the time-space-domain dispersion relation, the high-order FD coefficients are generated by using Taylor expansion and least squares. Numerical analyses and modeling examples demonstrate that our proposed schemes have higher accuracy and better stability than other conventional schemes, and thus larger time steps can be used to improve the computational efficiency in 3D case.

Compressibility in lattice Boltzmann on standard stencils: effects of deviation from reference temperature.

With growing interest in the simulation of compressible flows using the lattice Boltzmann (LB) method, a number of different approaches have been developed. These methods can be classified as pertaining to one of two major categories: (i) solvers relying on high-order stencils recovering the Navier-Stokes-Fourier equations, and (ii) approaches relying on classical first-neighbour stencils for the compressible Navier-Stokes equations coupled to an additional (LB-based or classical) solver for the energy balance equation. In most cases, the latter relies on a thermal Hermite expansion of the continuous equilibrium distribution function (EDF) to allow for compressibility. Even though recovering the correct equation of state at the Euler level, it has been observed that deviations of local flow temperature from the reference can result in instabilities and/or over-dissipation. The aim of the present study is to evaluate the stability domain of different EDFs, different collision models, with and without the correction terms for the third-order moments. The study is first based on a linear von Neumann analysis. The correction term for the space- and time-discretized equations is derived via a Chapman-Enskog analysis and further corroborated through spectral dispersion-dissipation curves. Finally, a number of numerical simulations are performed to illustrate the proposed theoretical study. This article is part of the theme issue 'Fluid dynamics, soft matter and complex systems: recent results and new methods'.

Accurate numerical solutions of 2-D elastodynamics problems using compact high-order stencils

A new approach for the increase in the order of accuracy of high-order numerical techniques used for 2-D time-dependent elasticity (structural dynamics and wave propagation) has been suggested on uniform square and rectangular meshes. It is based on the optimization of the coefficients of the corresponding semi-discrete stencil equations with respect to the local truncation error. It is shown that the second order of accuracy of the linear finite elements with 9-point stencils is optimal and cannot be improved. However, we have calculated new 9-point stencils with the second order of accuracy that yield more accurate results than those obtained by the linear finite elements. We have also developed new 25-point stencils (similar to those for the quadratic finite and isogeometric elements) with the optimal sixth order of accuracy for the non-diagonal mass matrix and with the optimal fourth order of accuracy for the diagonal mass matrix. The numerical results are in good agreement with the theoretical findings as well as they show a big increase in accuracy of the new stencils compared with those for the high-order finite elements. At the same number of degrees of freedom, the new approach yields significantly more accurate results than those obtained by the high-order (up to the eighth order) finite elements with the non-diagonal mass matrix. The numerical experiments also show that the new approach with 25-point stencils yields very accurate results for nearly incompressible materials with Poisson’s ratio 0.4999.

Compact high-order stencils with optimal accuracy for numerical solutions of 2-D time-independent elasticity equations

A new approach for the increase in the order of accuracy of high-order numerical techniques used for time-independent elasticity is suggested on uniform square and rectangular meshes. It is based on the optimization of the coefficients of the corresponding discrete stencil equations with respect to the local truncation error. It is shown that the second order of accuracy of 2-D linear finite elements with 9-point stencils is optimal and cannot be improved. However, the order of accuracy of 25-point stencils (similar to those for quadratic finite and isogeometric elements) can be significantly improved. We have developed new 25-point stencils for 2-D elastic problems with optimal 10th order of accuracy. The numerical results are in good agreement with the theoretical findings as well as they show a big increase in accuracy of the new stencils compared with those for high-order finite elements. At the same number of degrees of freedom, the new approach yields significantly more accurate results than those obtained by high-order (up to the tenth order) finite elements. The numerical experiments also show that the new approach with 25-point stencils yields very accurate results for nearly incompressible materials with Poisson’s ratio 0.4999.

Extreme-Scale High-Order WENO Simulations of 3-D Detonation Wave with 10 Million Cores

High-order stencil computations, frequently found in many applications, pose severe challenges to emerging many-core platforms due to the complexities of hardware architectures as well as the sophisticated computing and data movement patterns. In this article, we tackle the challenges of high-order WENO computations in extreme-scale simulations of 3D gaseous waves on Sunway TaihuLight. We design efficient parallelization algorithms and present effective optimization techniques to fully exploit various parallelisms with reduced memory footprints, enhanced data reuse, and balanced computation load. Test results show the optimized code can scale to 9.98 million cores, solving 12.74 trillion unknowns with 23.12 Pflops double-precision performance.

Optimization of Finite-Differencing Kernels for Numerical Relativity Applications

A simple optimization strategy for the computation of 3D finite-differencing kernels on many-cores architectures is proposed. The 3D finite-differencing computation is split direction-by-direction and exploits two level of parallelism: in-core vectorization and multi-threads shared-memory parallelization. The main application of this method is to accelerate the high-order stencil computations in numerical relativity codes. Our proposed method provides substantial speedup in computations involving tensor contractions and 3D stencil calculations on different processor microarchitectures, including Intel Knight Landing.

Mimetic discretization of the Eikonal equation with Soner boundary conditions

Motivated by a specific application in seismic reflection, the goal of this paper is to present a modified version of the Castillo–Grone mimetic gradient operators that allows for a high-order accurate solution of the Eikonal equation with Soner boundary conditions. The modified gradient operators utilize a non-staggered grid. In dimensions other than 1D, the modified gradient operators are expressed as Kronecker products of their corresponding 1D versions and some identity matrices. It is shown, that these modified 1D gradient operators are as accurate as the original gradient operators in terms of approximating first-order partial derivatives. It turns out, that in 1D one requires to solve two linear systems for finding a numerical solution of the Eikonal equation. Some examples show that the solution obtained by utilizing the modified operators increases its accuracy when incrementing the order of their approximation, something that does not occur when using the original operators. An iterative scheme is presented for the nonlinear 2D case. The method is of a quasi-Newton-like nature. At each iteration a linear system is built, with progressively higher-order stencils. The solution of the Fast Marching method is the initial guess. Numerical evidence indicates that high-order accurate solutions can be achieved.

Register optimizations for stencils on GPUs

The recent advent of compute-intensive GPU architecture has allowed application developers to explore high-order 3D stencils for better computational accuracy. A common optimization strategy for such stencils is to expose sufficient data reuse by means such as loop unrolling, with the expectation of register-level reuse. However, the resulting code is often highly constrained by register pressure. While current state-of-the-art register allocators are satisfactory for most applications, they are unable to effectively manage register pressure for such complex high-order stencils, resulting in sub-optimal code with a large number of register spills. In this paper, we develop a statement reordering framework that models stencil computations as a DAG of trees with shared leaves, and adapts an optimal scheduling algorithm for minimizing register usage for expression trees. The effectiveness of the approach is demonstrated through experimental results on a range of stencils extracted from application codes.

Methods for compressible fluid simulation on GPUs using high-order finite differences

We focus on implementing and optimizing a sixth-order finite-difference solver for simulating compressible fluids on a GPU using third-order Runge–Kutta integration. Since graphics processing units perform well in data-parallel tasks, this makes them an attractive platform for fluid simulation. However, high-order stencil computation is memory-intensive with respect to both main memory and the caches of the GPU. We present two approaches for simulating compressible fluids using 55-point and 19-point stencils. We seek to reduce the requirements for memory bandwidth and cache size in our methods by using cache blocking and decomposing a latency-bound kernel into several bandwidth-bound kernels. Our fastest implementation is bandwidth-bound and integrates 343 million grid points per second on a Tesla K40t GPU, achieving a 3.6× speedup over a comparable hydrodynamics solver benchmarked on two Intel Xeon E5-2690v3 processors. Our alternative GPU implementation is latency-bound and achieves the rate of 168 million updates per second.

Performance of the Double Absorbing Boundary Method when Applied to the 3D Acoustic Wave Equation

SUMMARY The double absorbing boundary (DAB) is a new highorder absorbing boundary condition for the scalar acoustic wave equation. It suppresses scattered waves at the edge of a boundary layer in computational domain boundary by using destructive interference analogous to a noise-cancelling headphone. This method has advantages in that it addresses some of the shortfalls in existing boundary conditions, such as the need for tuning in Perfectly Matched Layers or complex formulations at corners such as in high-order absorbing boundary conditions. We extend the original formulation of the DAB to three dimensions and higher-order stencils. Through numerical simulation we test the performance of the DAB by comparison with a reflecting boundary. We find that the DABC is a broadband attenuator with a power attenuation of 20-30dB using only six boundary cells. Increasing the order of the method improves accuracy for wavelengths less than 10 cells, whereas increasing the layer width does not improve accuracy. The method shows promise as a robust and computationally efficient boundary condition for seismic applications.

A framework for enhancing data reuse via associative reordering

The freedom to reorder computations involving associative operators has been widely recognized and exploited in designing parallel algorithms and to a more limited extent in optimizing compilers. In this paper, we develop a novel framework utilizing the associativity and commutativity of operations in regular loop computations to enhance register reuse. Stencils represent a particular class of important computations where the optimization framework can be applied to enhance performance. We show how stencil operations can be implemented to better exploit register reuse and reduce load/stores. We develop a multi-dimensional retiming formalism to characterize the space of valid implementations in conjunction with other program transformations. Experimental results demonstrate the effectiveness of the framework on a collection of high-order stencils.

Multigrid lattice Boltzmann method for accelerated solution of elliptic equations

A new solver for second-order elliptic partial differential equations (PDEs) based on the lattice Boltzmann method (LBM) and the multigrid (MG) technique is presented. Several benchmark elliptic equations are solved numerically with the inclusion of multiple grid-levels in two-dimensional domains at an optimal computational cost within the LB framework. The results are compared with the corresponding analytical solutions and numerical solutions obtained using the Stone's strongly implicit procedure. The classical PDEs considered in this article include the Laplace and Poisson equations with Dirichlet boundary conditions, with the latter involving both constant and variable coefficients. A detailed analysis of solution accuracy, convergence and computational efficiency of the proposed solver is given. It is observed that the use of a high-order stencil (for smoothing) improves convergence and accuracy for an equivalent number of smoothing sweeps. The effect of the type of scheduling cycle (V- or W-cycle) on the performance of the MG–LBM is analyzed. Next, a parallel algorithm for the MG–LBM solver is presented and then its parallel performance on a multi-core cluster is analyzed. Lastly, a practical example is provided wherein the proposed elliptic PDE solver is used to compute the electro-static potential encountered in an electro-chemical cell, which demonstrates the effectiveness of this new solver in complex coupled systems. Several orders of magnitude gains in convergence and parallel scaling for the canonical problems, and a factor of 5 reduction for the multiphysics problem are achieved using the MG–LBM.

HIGH-ORDER LATTICE-BOLTZMANN EQUATIONS AND STENCILS FOR MULTIPHASE MODELS

The lattice Boltzmann (LB) method, based on mesoscopic modeling of transport phenomena, appears to be an attractive alternative for the simulation of complex fluid flows. Examples of such complex dynamics are multiphase and multicomponent flows for which several LB models have already been proposed. However, due to theoretical or numerical reasons, some of these models may require application of high-order lattice-Boltzmann equations (LBEs) and stencils. Here, we will present a derivation of LBEs from the discrete-velocity Boltzmann equation (DVBE). By using the method of characteristics, high-order accurate equations are conveniently formulated with standard numerical methods for ordinary differential equations (ODEs). In particular, we will derive implicit LB schemes due to their stability properties. A simple algorithm is presented which enables implementation of the implicit schemes without resorting to, e.g. change of variables. Finally, some numerical experiments with high-order equations and stencils together with two specific multiphase models are presented.

Hierarchical parallelization and optimization of high-order stencil computations on multicore clusters

We present a scalable parallelization scheme for high-order stencil computations that also optimizes memory behavior on multicore clusters. Our multilevel approach combines: (i) inter-node parallelization via spatial decomposition; (ii) inter-core parallelization via multithreading and explicit non-uniform memory access (NUMA) control; (iii) data locality optimizations through auto-tuned tiling for efficient use of hierarchical memory; and (iv) register blocking and data parallelism via single-instruction multiple-data techniques to utilize registers and exploit data locality. The scheme is applied to a sixth-order stencil based finite-difference time-domain code. Weak-scaling parallel efficiency is over 98 % on 32,768 BlueGene/P processors. Multithreading with explicit NUMA control attains 9.9-fold speedup on a dual 12-core AMD Opteron system. Data locality optimizations achieve 7.7-fold reduction of the last level cache miss rate of Intel Nehalem, whereas register blocking increases data parallelism and thereby achieves 5.9 Gflops performance on a single core. Register blocking + multithreading optimizations achieve 5.8-fold speedup on a single quadcore Nehalem.

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.