High-order Finite Difference Discretization Research Articles

A unified Quasi-Spectral Viscosity (QSV) approach to shock capturing and large-eddy simulation

The Quasi-Spectral Viscosity (QSV) method is a novel closure for a high-order finite-difference discretization of the filtered compressible Navier-Stokes equations capable of unifying dynamic subfilter scale (SFS) modeling and shock capturing under a single mathematical framework. Its innovation lies in the introduction of a physical-space implementation of a spectral-like SFS dissipation term by leveraging residuals of filter operations, achieving two goals: (1) estimating the energy of the resolved solution near the grid cutoff; (2) imposing a plateau-cusp shape to the spectral distribution of the added dissipation. The QSV approach has been tested in a variety of flows to showcase its capability to act interchangeably as: a shock capturing method, in the Shu-Osher, shock/vortex or shock/wall interactions problems; or as a SFS closure, in subsonic Taylor Green Vortex (TGV), and supersonic/hypersonic turbulent channel flows. QSV performs well compared to previous eddy-viscosity closures and shock capturing methods in such test cases. In a supersonic TGV flow, a case which exhibits shock/turbulence interactions, QSV alone outperforms the simple superposition of separate numerical treatments for SFS turbulence and shocks. QSV's combined capability of simulating shocks and turbulence independently, as well as simultaneously, effectively achieves the unification of shock capturing and Large-Eddy Simulation.

Stability analysis of high-order finite-difference discretizations of the linearized forward-speed seakeeping problem

A high-order finite-difference method solution of the linearized, potential flow, seakeeping problem for a ship at steady forward speed was recently presented by Amini-Afshar et al. [1,2]. In this paper, we provide a detailed matrix-based eigenvalue stability analysis of this model, highlighting the sources of instability and the effects of possible remedies. In particular, we illustrate how both boundary treatment and grid stretching are important factors which are not typically captured by a von Neumann-type analysis. The new analysis shows that when grid stretching is used together with centered finite difference schemes, the method is generally unstable. The source of the instability can in some cases be traced to an effective downwinding of the convective terms. Stable solutions can be obtained either by introducing upwind-biased schemes for computing the convective derivatives on the free-surface, or by application of a mild filter at each time-step. A second source of instability is associated with the treatment of the convective derivatives of the free-surface elevation at points close to the domain boundaries. Here it is necessary to consider whether the surrounding fluid points lie in an upwind or a downwind direction. For upwinded points, ordinary one-sided differencing can be used, but for downwinded points we instead impose a Neumann-type boundary condition derived from the body and free-surface boundary conditions. As an example application to complement those already given in [1,2], the method is applied to solve the steady wave resistance problem and comparison is made to reference solutions for a two-dimensional floating cylinder and a submerged sphere. Estimates of the wave resistance of the Wigley hull are also compared with experimental measurements.

SQDFT: Spectral Quadrature method for large-scale parallel [formula omitted] Kohn–Sham calculations at high temperature

We present SQDFT: a large-scale parallel implementation of the Spectral Quadrature (SQ) method for O(N) Kohn–Sham Density Functional Theory (DFT) calculations at high temperature. Specifically, we develop an efficient and scalable finite-difference implementation of the infinite-cell Clenshaw–Curtis SQ approach, in which results for the infinite crystal are obtained by expressing quantities of interest as bilinear forms or sums of bilinear forms, that are then approximated by spatially localized Clenshaw–Curtis quadrature rules. We demonstrate the accuracy of SQDFT by showing systematic convergence of energies and atomic forces with respect to SQ parameters to reference diagonalization results, and convergence with discretization to established planewave results, for both metallic and insulating systems. We further demonstrate that SQDFT achieves excellent strong and weak parallel scaling on computer systems consisting of tens of thousands of processors, with near perfect O(N) scaling with system size and wall times as low as a few seconds per self-consistent field iteration. Finally, we verify the accuracy of SQDFT in large-scale quantum molecular dynamics simulations of aluminum at high temperature. Program summaryProgram Title: SQDFTProgram Files doi:http://dx.doi.org/10.17632/h4jnrmh9v3.1Licensing provisions: GNU General Public License 3 (GPL)Programming language: C/C++External routines/libraries: MVAPICH2 2.1 (http://mvapich.cse.ohio-state.edu/) or later.Nature of problem: Quantum mechanical calculation of static and dynamic properties of condensed matter at high temperature in the framework of Kohn–Sham Density Functional Theory (DFT).Solution method: High-order finite-difference discretization. Calculation of the electronic ground state using the self-consistent field (SCF) iteration in conjunction with Periodic Pulay extrapolation/mixing. Local reformulation of the electrostatics in terms of the electrostatic potential and pseudocharge densities. Solution of the Poisson equation using the Alternating Anderson–Richardson (AAR) method. Calculation of the electron density, band structure energy, electronic entropy energy, and atomic forces using the infinite-cell Clenshaw–Curtis Spectral Quadrature (SQ) approach, in which results for the infinite crystal are obtained by expressing quantities as bilinear forms or sums of bilinear forms, that are then approximated by spatially localized Clenshaw–Curtis quadrature rules. NVE molecular dynamics using the leapfrog method and NVT molecular dynamics using the Verlet algorithm with the Nose–Hoover thermostat. Extrapolation of electron density between molecular dynamics steps. Parallelization via domain decomposition.Restrictions: Cubical domain. Local Density Approximation (LDA). Troullier–Martins pseudopotentials without relativistic or non-linear core corrections.

SPARC: Accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory: Extended systems

As the second component of SPARC (Simulation Package for Ab-initio Real-space Calculations), we present an accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory (DFT) for extended systems. Specifically, employing a local formulation of the electrostatics, the Chebyshev polynomial filtered self-consistent field iteration, and a reformulation of the non-local force component, we develop a finite-difference framework wherein both the energy and atomic forces can be efficiently calculated to within desired accuracies in DFT. We demonstrate using a wide variety of materials systems that SPARC achieves high convergence rates in energy and forces with respect to spatial discretization to reference plane-wave result; exponential convergence in energies and forces with respect to vacuum size for slabs and wires; energies and forces that are consistent and display negligible ‘egg-box’ effect; accurate properties of crystals, slabs, and wires; and negligible drift in molecular dynamics simulations. We also demonstrate that the weak and strong scaling behavior of SPARC is similar to well-established and optimized plane-wave implementations for systems consisting up to thousands of electrons, but with a significantly reduced prefactor. Overall, SPARC represents an attractive alternative to plane-wave codes for performing DFT simulations of extended systems. Program summaryProgram title: SPARCCatalogue identifier: AFBR_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AFBR_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: GNU GPL v3No. of lines in distributed program, including test data, etc.: 93822No. of bytes in distributed program, including test data, etc.: 1386659Distribution format: tar.gzProgramming language: C/C++.Computer: Any system with C/C++ compiler.Operating system: Linux.RAM: Problem dependent. Ranges from 80 GB to 800 GB for a system with 2500 electrons.Classification: 7.3.External routines: PETSc 3.5.3 (http://www.mcs.anl.gov/petsc), MKL 11.2 (https://software.intel.com/en-us/intel-mkl), and MVAPICH2 2.1 (http://mvapich.cse.ohio-state.edu/).Does the new version supersede the previous version?: YesNature of problem: Calculation of the static and dynamic properties of isolated and extended systems in the framework of Kohn–Sham Density Functional Theory (DFT).Solution method:High-order finite-difference discretization. Local reformulation of the electrostatics in terms of the electrostatic potential and pseudocharge densities. Application of Bloch-periodic and zero-Dirichlet boundary conditions on the orbitals in the direction of periodicity and vacuum, respectively. Application of periodic and Dirichlet boundary conditions on the electrostatic potential in the direction of periodicity and vacuum, respectively. Integration over the Brillouin zone for extended systems using the Monkhorst–Pack grid. Calculation of the electronic ground-state using the Chebyshev polynomial filtered self-consistent field iteration in conjunction with Anderson based extrapolation/mixing schemes. Reformulation of the non-local component of the force. Geometry optimization using the Polak–Ribiere variant of non-linear conjugate gradients with secant line search. NVE molecular dynamics using the leapfrog method. Parallelization via domain decomposition and over Brillouin zone integration.Reasons for new version:To enable the study of extended systems like crystals, slabs, and wires using SPARC.Summary of revisions:Incorporated the ability to study the static and dynamic properties of crystals, slabs, and wires.Restrictions:System size less than ∼4000 electrons. Local Density Approximation (LDA). Troullier–Martins pseudopotentials without relativistic or non-linear core corrections. Domain has to be cuboidal.Running time:Problem dependent. Timing results for selected examples provided in the paper.

SPARC: Accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory: Isolated clusters

As the first component of SPARC (Simulation Package for Ab-initio Real-space Calculations), we present an accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory (DFT) for isolated clusters. Specifically, utilizing a local reformulation of the electrostatics, the Chebyshev polynomial filtered self-consistent field iteration, and a reformulation of the non-local component of the force, we develop a framework using the finite-difference representation that enables the efficient evaluation of energies and atomic forces to within the desired accuracies in DFT. Through selected examples consisting of a variety of elements, we demonstrate that SPARC obtains exponential convergence in energy and forces with domain size; systematic convergence in the energy and forces with mesh-size to reference plane-wave result at comparably high rates; forces that are consistent with the energy, both free from any noticeable ‘egg-box’ effect; and accurate ground-state properties including equilibrium geometries and vibrational spectra. In addition, for systems consisting up to thousands of electrons, SPARC displays weak and strong parallel scaling behavior that is similar to well-established and optimized plane-wave implementations, but with a significantly reduced prefactor. Overall, SPARC represents an attractive alternative to plane-wave codes for practical DFT simulations of isolated clusters.Program summaryProgram title: SPARCCatalogue identifier: AFBL_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AFBL_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: GNU GPL v3No. of lines in distributed program, including test data, etc.: 47525No. of bytes in distributed program, including test data, etc.: 826436Distribution format: tar.gzProgramming language: C/C++.Computer: Any system with C/C++ compiler.Operating system: Linux.RAM: Problem dependent. Ranges from 80 GB to 800 GB for a system with 2500 electrons.Classification: 7.3.External routines: PETSc 3.5.3 (http://www.mcs.anl.gov/petsc), MKL 11.2 (https://software.intel.com/en-us/intel-mkl), and MVAPICH2 2.1 (http://mvapich.cse.ohio-state.edu/).Nature of problem:Calculation of the electronic and structural ground-states for isolated clusters in the framework of Kohn–Sham Density Functional Theory (DFT).Solution method:High-order finite-difference discretization. Local reformulation of the electrostatics in terms of the electrostatic potential and pseudocharge densities. Calculation of the electronic ground-state using the Chebyshev polynomial filtered Self-Consistent Field (SCF) iteration in conjunction with Anderson extrapolation/mixing. Evaluation of boundary conditions for the electrostatic potential through a truncated multipole expansion. Reformulation of the non-local component of the force. Geometry optimization using the Polak–Ribiere variant of non-linear conjugate gradients with secant line search.Restrictions:System size less than ∼4000 electrons. Local Density Approximation (LDA). Troullier–Martins pseudopotentials without relativistic or non-linear core corrections.Running time:Problem dependent. Timing results for selected examples provided in the paper.

High Order Finite Difference Discretization for Composite Grid Hierarchy and Its Applications

AbstractWe introduce efficient approaches to construct high order finite difference discretizations for solving partial differential equations, based on a composite grid hierarchy. We introduce a modification of the traditional point clustering algorithm, obtained by adding restrictive parameters that control the minimal patch length and the size of the buffer zone. As a result, a reduction in the number of interfacial cells is observed. Based on a reasonable geometric grid setting, we discuss a general approach for the construction of stencils in a composite grid environment. The straightforward approach leads to an ill-posed problem. In our approach we regularize this problem, and transform it into solving a symmetric system of linear of equations. Finally, a stencil repository has been designed to further reduce computational overhead. The effectiveness of the discretizations is illustrated by numerical experiments on second order elliptic differential equations.

Efficiency Benchmarking of an Energy Stable High-Order Finite Difference Discretization

In this paper, results are presented for a number of benchmark cases, proposed at the 2nd International Workshop on High-Order CFD Methods in Cologne, Germany, in 2013. A robust high-order-accurate finite difference method was used that was developed during the last 10–15 years. The robustness stems from the fact that the entire numerical procedure, including various kinds of boundary conditions, can be proven stable. This paper outlines the methodology, and it summarizes results presented at the workshop along with some more data and tests.

Photorealistic visualization and fluid animation: coupling of Maya with a two-phase Navier-Stokes fluid solver

We have coupled the three-dimensional solver for the two-phase incompressible Navier-Stokes equations NaSt3DGPF with Autodesk Maya. Maya is the industry standard software framework for the creation of three-dimensional animations. The parallel level-set-based fluid solver NaSt3DGPF simulates the interaction of two fluids like air and water. It uses high-order finite difference discretization methods that are designed for physics applications. By coupling both applications, we are now able to set up scientific fluid simulations in an easy-to-use user interface. Moreover, the rendering techniques provided by Maya allow us to create photorealistic visualizations for computational fluid dynamics problems and support the creation of highly visually realistic fluid simulations for animation movies. Altogether, we obtain an easy usable and fully coupled fluid animation toolkit for two-phase fluid simulations. These are the first published results of the full integration of a physics-oriented, high-order grid-based parallel two-phase fluid solver in Maya, at least to our knowledge.

Stable and Accurate Interpolation Operators for High-Order Multiblock Finite Difference Methods

Block-to-block interface interpolation operators are constructed for several common high-order finite difference discretizations. In contrast to conventional interpolation operators, these new interpolation operators maintain the strict stability, accuracy, and conservation properties of the base scheme even when nonconforming grids or dissimilar operators are used in adjoining blocks. The stability properties of the new operators are verified using eigenvalue analysis, and the accuracy properties are verified using numerical simulations of the Euler equations in two spatial dimensions.

Finite difference methods for second order in space, first order in time hyperbolic systems and the linear shifted wave equation as a model problem in numerical relativity

Motivated by the problem of solving the Einstein equations, we discuss high order finite difference discretizations of first order in time, second order in space hyperbolic systems.Particular attention is paid to the case when first order derivatives that can be identified with advection terms are approximated with non-centered finite difference operators.We first derive general properties of these discrete operators, then we extend a known result on numerical stability for such systems to general order of accuracy.As an application we analyze the shifted wave equation, including the behavior of the numerical phase and group speeds at different orders of approximations. Special attention is paid to when the use of off-centered schemes improves the accuracy over the centered schemes.

Stable Boundary Treatment for the Wave Equation on Second-Order Form

A stable and accurate boundary treatment is derived for the second-order wave equation. The domain is discretized using narrow-diagonal summation by parts operators and the boundary conditions are imposed using a penalty method, leading to fully explicit time integration. This discretization yields a stable and efficient scheme. The analysis is verified by numerical simulations in one-dimension using high-order finite difference discretizations, and in three-dimensions using an unstructured finite volume discretization.

Stable and accurate wave-propagation in discontinuous media

A time stable discretization is derived for the second-order wave equation with discontinuous coefficients. The discontinuity corresponds to inhomogeneity in the underlying medium and is treated by splitting the domain. Each (homogeneous) sub domain is discretized using narrow-diagonal summation by parts operators and, then, patched to its neighbors by using a penalty method, leading to fully explicit time integration. This discretization yields a time stable and efficient scheme. The analysis is verified by numerical simulations in one-dimension using high-order finite difference discretizations, and in three-dimensions using an unstructured finite volume discretization.

High-order local absorbing conditions for the wave equation: Extensions and improvements

The solution of the time-dependent wave equation in an unbounded domain is considered. An artificial boundary B is introduced which encloses a finite computational domain. On B an absorbing boundary condition (ABC) is imposed. A formulation of local high-order ABCs recently proposed by Hagstrom and Warburton and based on a modification of the Higdon ABCs, is further developed and extended in a number of ways. First, the ABC is analyzed in new ways and important information is extracted from this analysis. Second, The ABCs are extended to the case of a dispersive medium, for which the Klein–Gordon wave equation governs. Third, the case of a stratified medium is considered and the way to apply the ABCs to this case is explained. Fourth, the ABCs are extended to take into account evanescent modes in the exact solution. The analysis is applied throughout this paper to two-dimensional wave guides. Two numerical algorithms incorporating these ABCs are considered: a standard semi-discrete finite element formulation in space followed by time-stepping, and a high-order finite difference discretization in space and time. Numerical examples are provided to demonstrate the performance of the extended ABCs using these two methods.

Comparison of some upwind-biased high-order formulations with a second-order central-difference scheme for time integration of the incompressible Navier-Stokes equations

The paper compares high-order finite-difference discretizations with a second-order central-difference scheme for the time integration of the incompressible Navier-Stokes equations. Conservative and non-conservative forms of the convection terms are discretized using fifth-order accurate upwind-biased based approximations. The high-order non-conservative treatment of the convection terms exhibits the best accuracy at high resolutions but deteriorates rapidly as the resolution decreases. The high-order conservative treatment of the convection terms, despite the second-order finite-volume operator, exhibits high-order accuracy but does not show any advantage over the non-conservative formulation. Combinations of the divergence and gradient operators based on second and fourth-order accurate central-difference approximations are used to construct high-order Laplacians in the pressure equation. There is little evidence that the high-order treatment of the pressure equation adds sufficiently to the overall accuracy of the scheme to justify the extra computational effort associated with it. Simulations of turbulent channel flow indicate that the second-order central-difference scheme resolves the turbulent spectrum better than the high-order upwind schemes.

Low-Dissipation and Low-Dispersion Runge–Kutta Schemes for Computational Acoustics

In this paper, we investigate accurate and efficient time advancing methods for computational acoustics, where nondissipative and nondispersive properties are of critical importance. Our analysis pertains to the application of Runge–Kutta methods to high-order finite difference discretization. In many CFD applications, multistage Runge–Kutta schemes have often been favored for their low storage requirements and relatively large stability limits. For computing acoustic waves, however, the stability consideration alone is not sufficient, since the Runge–Kutta schemes entail both dissipation and dispersion errors. The time step is now limited by the tolerable dissipation and dispersion errors in the computation. In the present paper, it is shown that if the traditional Runge–Kutta schemes are used for time advancing in acoustic problems, time steps greatly smaller than those allowed by the stability limit are necessary. Low-dissipation and low-dispersion Runge–Kutta (LDDRK) schemes are proposed, based on an optimization that minimizes the dissipation and dispersion errors for wave propagation. Optimizations of both single-step and two-step alternating schemes are considered. The proposed LDDRK schemes are remarkably more efficient than the classical Runge–Kutta schemes for acoustic computations. Moreover, low storage implementations of the optimized schemes are discussed. Special issues of implementing numerical boundary conditions in the LDDRK schemes are also addressed.

Families of High Order Accurate Discretizations of Some Elliptic Problems

In this paper we construct and analyze high order finite difference discretizations of a class of elliptic partial differential equations. In particular, two one-parameter families of fourth order HODIE discretizations of the Helmholtz equation are derived and a discretization optimal with respect to a certain norm of the truncation error is identified. The use of compact nine-point formulas of positive type admits both fast direct methods and standard iterative methods for the solution of the resulting systems of linear equations. Extensions yielding sixth order accuracy for the Helmholtz equation and fourth order accuracy for a more general operator are given. Finally, numerical results demonstrating the effectiveness of the discretizations for a wide range of problems are presented.