Grid Stencil Research Articles

Grid-Characteristic Numerical Method on an Irregular Grid with Extending the Interpolation Stencil

A grid-characteristic numerical method for solving a multidimensional transport equation on an unstructured grid with an order higher than one is proposed; this method does not use auxiliary points on edges and faces. The avoidance of auxiliary points on edges and faces simplifies the topology of the computational grid during its motion, which is important for solving dynamic problems of mechanics of deformable solids. To increase the approximation order, an analog of the grid stencil extension implemented for an unstructured grid is used. Results of testing the proposed numerical scheme for continuously differentiable, continuous, discontinuous solutions are presented.

Modified Explicit Decoupled Group Scheme οn Helmholtz Equation

In this paper, the formulation of a new group iterative method called the Modified Explicit Decoupled Group method in solving the two dimensional Helmholtz equation is described. The method is derived using a combination of the five-point finite difference approximation on the rotated grid stencil together with the five-point centred difference approximation on the standard grid stencils. Numerical experimentations of this new formulation shows significant improvement in computational complexity and execution timings over the original Explicit Decoupled Group method [2].

Acceleration of high-order combined compact finite difference scheme for simulating three-dimensional flow and heat transfer problems in GPUs

In this article, the high-order upwinding combined compact difference scheme developed in a three-point grid stencil is applied to solve the incompressible Navier-Stokes (NS) and energy equations in three dimensions. The time integrator with symplectic property is employed to approximate the temporal derivative term in inviscid Euler equation so as to numerically retain the embedded Hamiltions and Casimir to get long-time accurate solutions. For the sake of computational efficiency in solving the three-dimensional NS equations, all the calculations will be accelerated using the hybrid CUDA and OpenAcc GPU programing models. The parallel speedup performance compared to the multicore of an Intel Xeon E5-2690V5 CPU is reported.

Accuracy analysis of gradient reconstruction on isotropic unstructured meshes and its effects on inviscid flow simulation

The accuracy of gradient reconstruction methods on unstructured meshes is analyzed both mathematically and numerically. Mathematical derivations reveal that, for gradient reconstruction based on the Green-Gauss theorem (the GG methods), if the summation of first-and-lower-order terms does not counterbalance in the discretized integral process, which rarely occurs, second-order accurate approximation of face midpoint value is necessary to produce at least first-order accurate gradient. However, gradient reconstruction based on the least-squares approach (the LSQ methods) is at least first-order on arbitrary unstructured grids. Verifications are performed on typical isotropic grid stencils by analyzing the relationship between the discretization error of gradient reconstruction and the discretization error of the face midpoint value approximation of a given analytic function. Meanwhile, the numerical accuracy of gradient reconstruction methods is examined with grid convergence study on typical isotropic grids. Results verify the phenomenon of accuracy degradation for the GG methods when the face midpoint value condition is not satisfied. The LSQ methods are proved to be at least first-order on all tested isotropic grids. To study gradient accuracy effects on inviscid flow simulation, solution errors are quantified using the Method of Manufactured Solutions (MMS) which was validated before adoption by comparing with an exact solution case, i.e., the 2-dimensional (2D) inviscid isentropic vortex. Numerical results demonstrate that the order of accuracy (OOA) of gradient reconstruction is crucial in determining the OOA of numerical solutions. Solution accuracy deteriorates seriously if gradient reconstruction does not reach first-order.

New very high-order upwind multi-layer compact (MLC) schemes with spectral-like resolution for flow simulations

Numerical simulations of multi-scale flow problems such as hypersonic boundary layer transition, turbulent flows, computational aeroacoustics and various other flow problems with complex physics require high-order methods with high spectral resolutions. For instance, the receptivity mechanisms in the hypersonic boundary layer are the resonant interactions between forcing waves and boundary-layer waves, and the complex wave interactions are difficult to be accurately predicted by conventional low-order numerical methods. High-order methods, which are robust and accurate in resolving a wide range of time and length scales, are required. Currently, the high-order finite difference methods for simulations of hypersonic flows are usually upwind schemes or compact schemes with fifth-order accuracy or lower [1]. The objective of this paper is to develop and analyze a new very high-order numerical scheme with the spectral-like resolution for flow simulations on structured grids, with focus on smooth flow problems involving multiple scales. Specifically, a new upwind multi-layer compact (MLC) scheme with spectral-like resolution up to seventh order is derived in a finite difference framework. By using the ‘multi-layer’ idea, which introduces first derivatives into the MLC schemes and approximates the second derivatives, the resolution of the MLC schemes can be significantly improved within a compact grid stencil. The auxiliary equations are required and they are the only nontrivial equations, which contributes to good computational efficiency. In addition, the upwind MLC schemes are derived based on the idea of constructing upwind schemes on centered stencils with adjustable parameters to control the dissipation. Fourier analysis is performed to show that the MLC schemes have small dissipation and dispersion in a very wide range of wavenumbers in both one- and two-dimensional cases, and the anisotropic error is much smaller than conventional finite difference methods in the two-dimensional case. Furthermore, the stability analysis with matrix method shows that high-order boundary closure schemes are stable because of compactness of the stencils. The accuracies and rates of convergence of the new schemes are validated by numerical experiments of the linear advection equation, the nonlinear Euler equations, and the Navier–Stokes equations. The numerical results show that good computational efficiency, very high-order accuracies, and high spectral resolutions especially on coarse meshes can be attained with the MLC schemes. Overall, the MLC scheme has the properties of simple formulations, high-order accuracies, spectral-like resolutions, and compact stencils, and it is suitable for accurate simulation of smooth multi-scale flows with complex physics.

Numerical modeling of lock-exchange gravity/turbidity currents by a high-order upwinding combined compact difference scheme

This study presents two-dimensional direct numerical simulations for sediment-laden current with higher density propagating forward through a lighter ambient water. The incompressible NavierStokes equations including the buoyancy force for the density difference between the light and heavy fluids are solved by a finite difference scheme based on a structured mesh. The concentration transport equations are used to explore such rich transport phenomena as gravity and turbidity currents. Within the framework of an Upwinding Combined Compact finite Difference (UCCD) scheme, rigorous determination of weighting coefficients underlies the modified equation analysis and the minimization of the numerical modified wavenumber. This sixth-order UCCD scheme is implemented in a four-point grid stencil to approximate advection and diffusion terms in the concentration transport equations and the first-order derivative terms in the Navier-Stokes equations, which can greatly enhance convective stability and increase dispersive accuracy at the same time. The initial discontinuous concentration field is smoothed by solving a newly proposed Heaviside function to prevent numerical instabilities and unreasonable concentration values. A two-step projection method is then applied to obtain the velocity field. The numerical algorithm shows a satisfying ability to capture the generation, development, and dissipation of the Kelvin-Helmholz instabilities and turbulent billows at the interface between the current and the ambient fluid. The simulation results also are compared with the data in published literatures and good agreements are found to prove that the present numerical model can well reproduce the propagation, particle deposition, and mixing processes of lock-exchange gravity and turbidity currents.

Comparison of dealiasing schemes in large-eddy simulation of neutrally stratified atmospheric flows

Abstract. Aliasing errors arise in the multiplication of partial sums, such as those encountered when numerically solving the Navier–Stokes equations, and can be detrimental to the accuracy of a numerical solution. In this work, a performance and cost analysis is proposed for widely used dealiasing schemes in large-eddy simulation, focusing on a neutrally stratified, pressure-driven atmospheric boundary-layer flow. Specifically, the exact 3∕2 rule, the Fourier truncation method, and a high-order Fourier smoothing method are intercompared. Tests are performed within a newly developed mixed pseudo-spectral finite differences large-eddy simulation code, parallelized using a two-dimensional pencil decomposition. A series of simulations are performed at varying resolution, and key flow statistics are intercompared among the considered runs and dealiasing schemes. The three dealiasing methods compare well in terms of first- and second-order statistics for the considered cases, with modest local departures that decrease as the grid stencil is reduced. Computed velocity spectra using the 3∕2 rule and the FS method are in good agreement, whereas the FT method yields a spurious energy redistribution across wavenumbers that compromises both the energy-containing and inertial sublayer trends. The main advantage of the FS and FT methods when compared to the 3∕2 rule is a notable reduction in computational cost, with larger savings as the resolution is increased (15 % for a resolution of 1283, up to a theoretical 30 % for a resolution of 20483).

A combined compact finite difference scheme for predicting the evolution of a mean curvature driven interface

This study is aimed to develop a high-order finite difference scheme to simulate the mean curvature equation for accurately predicting the time-evolving mean curvature driven interface. Within the semi-discretization framework, the optimal third-order accurate temporal scheme is applied for the approximation of the time derivative term. In a three-point grid stencil, a combined compact difference scheme has been developed for offering a fifth-order spatial accuracy for the first-order and second-order derivative terms shown in the level set equation for simulating mean curvature flow. In the simulation of the mean curvature equation, reinitialization procedure has been performed to improve the prediction of curve motion. The aim of this study is to give insights into the issue about “how a mean curvature driven curve is varied with time”. Two mechanisms leading to the change of curve slope are attributed partly to the damping Laplacian operator and partly to the embedded nonlinear differential operator. These mechanisms have been studied in detail with respect to the curvature of curve. The effect of performing reinitialization is also studied numerically.

Reprint of: Accuracy consideration by DRP schemes for DNS and LES of compressible flow computations

Several dispersion relation-preserving (DRP) spatially central discretizations are considered as the base scheme in the framework of the Yee & Sjögreen low dissipative nonlinear filter approach. In addition, the nonlinear filter of Yee & Sjögreen with shock-capturing and long time integration capabilities is used to replace the standard DRP linear filter for both smooth flows and flows containing discontinuities. DRP schemes for computational aeroacoustics (CAA) focus on dispersion error consideration for long time linear wave propagation rather than the formal order of accuracy of the scheme. The resulting DRP schemes usually have wider grid stencils and increased CPU operations count compared with standard central schemes of the same formal order of accuracy. For discontinuous initial data and long time wave propagation of smooth acoustic waves, various space and time DRP linear filter are needed. For acoustic waves interacting with shocks and turbulence induced noise, DRP schemes with linear filters alone usually are not capable of simulating such flows. The investigation presented in this paper is focused on the possible gain in efficiency and accuracy by spatial DRP schemes over standard central schemes having the same grid stencil width for general direct numerical simulations (DNS) and large eddy simulations (LES) of compressible flows. Representative test cases for both smooth flows and problems containing discontinuities for 3D DNS of compressible gas dynamics are included.

Numerical solutions to a two-component Camassa–Holm equation

Abstract In the present paper, a three-step iterative algorithm for solving a two-component Camassa–Holm (2CH) equation is presented. In the first step, the time-dependent equation for the horizontal fluid velocity with nonlinear convection is solved. Then an inhomogeneous Helmholtz equation is solved. Finally, the equation for modeling the transport of density is solved in the third step. The differential order of 2CH equation has been reduced in order to facilitate numerical scheme development in a comparatively smaller grid stencil. In this study, a fifth-order spatially accurate upwinding combined compact difference scheme (UCCD5) which differs from that in Sheu et al. (2011) is developed in a four-point grid stencil for approximating the first-order derivative term. For the purpose of retaining long-time Hamiltonians in the 2CH equation, the time integrator (or time-stepping scheme) chosen is symplectic. Various numerical experiments such as the single peakon, peakon–antipeakon interaction and dam-break problems are conducted to illustrate the effectiveness of the proposed numerical method. It is shown that both the Hamiltonians and Casimir functions are conserved well for all problems.

Accuracy consideration by DRP schemes for DNS and LES of compressible flow computations

Several dispersion relation-preserving (DRP) spatially central discretizations are considered as the base scheme in the framework of the Yee & Sjögreen low dissipative nonlinear filter approach. In addition, the nonlinear filter of Yee & Sjögreen with shock-capturing and long time integration capabilities is used to replace the standard DRP linear filter for both smooth flows and flows containing discontinuities. DRP schemes for computational aeroacoustics (CAA) focus on dispersion error consideration for long time linear wave propagation rather than the formal order of accuracy of the scheme. The resulting DRP schemes usually have wider grid stencils and increased CPU operations count compared with standard central schemes of the same formal order of accuracy. For discontinuous initial data and long time wave propagation of smooth acoustic waves, various space and time DRP linear filter are needed. For acoustic waves interacting with shocks and turbulence induced noise, DRP schemes with linear filters alone usually are not capable of simulating such flows. The investigation presented in this paper is focused on the possible gain in efficiency and accuracy by spatial DRP schemes over standard central schemes having the same grid stencil width for general direct numerical simulations (DNS) and large eddy simulations (LES) of compressible flows. Representative test cases for both smooth flows and problems containing discontinuities for 3D DNS of compressible gas dynamics are included.

A two-scale approach for efficient on-the-fly operator assembly in massively parallel high performance multigrid codes

Large scale matrix-free finite element implementations save memory and are often significantly faster than implementations using classical sparse matrix techniques. They are especially well suited for massively parallel geometric multigrid solvers in combination with hierarchical hybrid grids on polyhedral domains. In the case of constant coefficients, the number of different stencil entries depends only on the coarse grid and does not increase with the number of refinement levels. However, for non-polyhedral domains the situation changes. Then even for the Laplace operator, the element mapping leads to fine grid stencils that can vary from grid point to grid point. Traditional matrix-free techniques that are based on an element-wise assembly then result in a considerably increase in computational cost. To compensate for this shortcoming, we introduce a new two-scale approach that uses a surrogate operator. It exploits a piecewise polynomial approximation of the entries of the stencil of the fine grid operator with respect to the coarse mesh size. The low-cost evaluation of these surrogate polynomials results in an efficient stencil assembly on-the-fly for non-polyhedral domains. We discuss and illustrate numerically two-scale a priori bounds. The accuracy of the approximate solution can be further improved if combined with a double discretization technique. A careful performance analysis in combination with a hardware–aware code optimization based on the Execution–Cache–Memory model yields a significant speed up. Weak and strong scaling results illustrate the potential of this new two-scale approach within large scale PDE simulations.

Acceleration of Linear Finite-Difference Poisson-Boltzmann Methods on Graphics Processing Units.

Electrostatic interactions play crucial roles in biophysical processes such as protein folding and molecular recognition. Poisson-Boltzmann equation (PBE)-based models have emerged as widely used in modeling these important processes. Though great efforts have been put into developing efficient PBE numerical models, challenges still remain due to the high dimensionality of typical biomolecular systems. In this study, we implemented and analyzed commonly used linear PBE solvers for the ever-improving graphics processing units (GPU) for biomolecular simulations, including both standard and preconditioned conjugate gradient (CG) solvers with several alternative preconditioners. Our implementation utilizes the standard Nvidia CUDA libraries cuSPARSE, cuBLAS, and CUSP. Extensive tests show that good numerical accuracy can be achieved given that the single precision is often used for numerical applications on GPU platforms. The optimal GPU performance was observed with the Jacobi-preconditioned CG solver, with a significant speedup over standard CG solver on CPU in our diversified test cases. Our analysis further shows that different matrix storage formats also considerably affect the efficiency of different linear PBE solvers on GPU, with the diagonal format best suited for our standard finite-difference linear systems. Further efficiency may be possible with matrix-free operations and integrated grid stencil setup specifically tailored for the banded matrices in PBE-specific linear systems.

A high-order finite volume method for solving one-dimensional convection and diffusion equations

ABSTRACTA spatially high-order finite volume method for solving convection and diffusion equations is developed and tested in this work. The method performs a high-order piecewise polynomial reconstruction of the local flow field based on the relationship between Taylor’s series expansion and the volume-averaged flow quantities. A 5 × 5 matrix inversion for each cell is done to compute the cell-center variables and derivatives up to fourth order. While a fixed symmetric grid stencil is maintained in smooth flow regions, a detector for large change in linear data slopes is developed to trigger the use of ENO stencil around flow discontinuities. Regular time integration scheme such as the four-stage Runge–Kutta method or the Euler implicit method is used for time integration. The present finite volume method is shown to be spatially fifth-order accurate for the linear convection equation, fourth-order accurate for the linear diffusion equation, and fourth-order accurate for the linear convection–diffusion equation. The shocks captured in solving the inviscid Burger’s equation are sharp and oscillation free. For the system of Euler equations, a characteristic limiter is further developed to limit the growth of total variation of the solution. Test examples solving shock-tube problems and the interactions of two blast waves show that various flow discontinuities are captured sharply without spurious oscillations.

Simulation of Maxwell's Equations on GPU Using a High-Order Error-Minimized Scheme

AbstractIn this study an explicit Finite Difference Method (FDM) based scheme is developed to solve the Maxwell's equations in time domain for a lossless medium. This manuscript focuses on two unique aspects – the three dimensional time-accurate discretization of the hyperbolic system of Maxwell equations in three-point non-staggered grid stencil and it's application to parallel computing through the use of Graphics Processing Units (GPU). The proposed temporal scheme is symplectic, thus permitting conservation of all Hamiltonians in the Maxwell equation. Moreover, to enable accurate predictions over large time frames, a phase velocity preserving scheme is developed for treatment of the spatial derivative terms. As a result, the chosen time increment and grid spacing can be optimally coupled. An additional theoretical investigation into this pairing is also shown. Finally, the application of the proposed scheme to parallel computing using one Nvidia K20 Tesla GPU card is demonstrated. For the benchmarks performed, the parallel speedup when compared to a single core of an Intel i7-4820K CPU is approximately 190x.

Development of a Combined Compact Difference Scheme to Simulate Soliton Collision in a Shallow Water Equation

AbstractIn this paper a three-step scheme is applied to solve the Camassa-Holm (CH) shallow water equation. The differential order of the CH equation has been reduced in order to facilitate development of numerical scheme in a comparatively smaller grid stencil. Here a three-point seventh-order spatially accurate upwinding combined compact difference (CCD) scheme is proposed to approximate the firstorder derivative term. We conduct modified equation analysis on the CCD scheme and eliminate the leading discretization error terms for accurately predicting unidirectional wave propagation. The Fourier analysis is carried out as well on the proposed numerical scheme to minimize the dispersive error. For preserving Hamiltonians in Camassa- Holm equation, a symplecticity conserving time integrator has been employed. The other main emphasis of the present study is the use of u–P–α formulation to get nondissipative CH solution for peakon-antipeakon and soliton-anticuspon head-on wave collision problems.

Simulation of Incompressible Free Surface Flow Using the Volume Preserving Level Set Method

AbstractThis study aims to develop a numerical scheme in collocated Cartesian grids to solve the level set equation together with the incompressible two-phase flow equations. A seventh-order accurate upwinding combined compact difference (UCCD7) scheme has been developed for the approximation of the first-order spatial derivative terms shown in the level set equation. Developed scheme has a higher accuracy with a three-point grid stencil to minimize phase error. To preserve the mass of each phase all the time, the temporal derivative term in the level set equation is approximated by the sixth-order accurate symplectic Runge-Kutta (SRK6) scheme. All the simulated results for the dam-break, Rayleigh-Taylor instability, bubble rising, two-bubble merging, and milkcrown problems in two and three dimensions agree well with the available numerical or experimental results.

New explicit group iterative methods in the solution of three dimensional hyperbolic telegraph equations

In this paper, new group iterative numerical schemes based on the centred and rotated (skewed) seven-point finite difference discretisations are proposed for the solution of a three dimensional second order hyperbolic telegraph equation, subject to specific initial and Dirichlet boundary conditions. Both schemes are shown to be of second order accuracies and unconditionally stable. The scheme derived from the rotated grid stencil results in a reduced linear system with lower computational complexity compared to the scheme derived from the centred approximation formula. A comparative study with other common point iterative methods based on the seven-point centred difference approximation together with their computational complexity analyses is also presented.

A higher order perfectly matched layer formulation for finite-difference time-domain seismic wave modeling

We have developed a higher order perfectly matched layer (PML) formulation to improve the absorption performance for finite-difference time-domain seismic modeling. First, we outlined a new unsplit “correction” approach, which allowed for traditional, first-order PMLs to be added directly to existing codes in a straightforward manner. Then, using this framework, we constructed a PML formulation that can be used to construct higher order PMLs of arbitrary order. The greater number of degrees of freedom associated with the higher order PML allow for enhanced flexibility of the PML stretching functions, thus potentially facilitating enhanced absorption performance. We found that the new approach can offer increased elastodynamic absorption, particularly for evanescent waves. We also discovered that the extra degrees of freedom associated with the higher order PML required careful optimization if enhanced absorption was to be achieved. Furthermore, these extra degrees of freedom increased the computational requirements in comparison with first-order schemes. We reached our formulations using one compact equation thus increasing the ease of implementation. Additionally, the formulations are based on a recursive integration approach that reduce PML memory requirements, and do not require special consideration for corner regions. We tested the new formulations to determine their ability to absorb body waves and surface waves. We also tested standard staggered grid stencils and rotated staggered grid stencils.

Numerical investigation of chemotaxic phenomenon in incompressible viscous fluid flow

This study aims to get some details of the chemotaxic phenomena in a domain containing an incompressible viscous fluid by the proposed numerical method. We are also aimed to numerically revisit the classical blow-up phenomenon regarding the initially prescribed cell density in an infinitely larger domain. To get a computationally more accurate scheme for solving the Keller-Segel (KS) [1] equations with/without coupling with the nonlinear hydrodynamic equations, a combined compact difference scheme of fifth-order spatial accuracy is developed in a three-point grid stencil. Three different sets of Keller-Segel equations, which are all amenable to exact solutions with/without coupling with the incompressible Navier–Stokes (NS) equations, are solved to numerically verify the proposed combined compact difference scheme and the incompressible flow solver. How the cell density and chemical concentration are affected by flow convection and diffusion effects will be investigated by solving the coupled KS–NS differential equations.