Stencil Size Research Articles

Prefactored optimized compact finite-difference schemes for second spatial derivatives

We have described systematically the processes of developing prefactored optimized compact schemes for second spatial derivatives. First, instead of emphasizing high resolution of a single monochromatic wave, we focus on improving the representation of the compact finite difference schemes over a wide range of wavenumbers. This leads to the development of the optimized compact schemes whose coefficients will be determined by Fourier analysis and the least-squares optimization in the wavenumber domain. The resulted optimized compact schemes provide the maximum resolution in spatial directions for the simulation of wave propagations. However, solving for each spatial derivative using these compact schemes requires the inversion of a band matrix. To resolve this issue, we propose a prefactorization strategy that decomposes the original optimized compact scheme into forward and backward biased schemes, which can be solved explicitly. We achieve this by ensuring a property that the real numerical wavenumbers of both the forward and backward biased stencils are the same as that of the original central compact scheme, and their imaginary numerical wavenumbers have the same values but with opposite signs. This property guarantees that the original optimized compact scheme can be completely recovered after the summation of the forward and backward finite difference operators. These prefactored optimized compact schemes have smaller stencil sizes than even those of the original compact schemes, and hence, they can take full advantage of the computer caches without sacrificing their resolving power. Comparisons were made throughout with other well-known schemes.

An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Convex Hamilton-Jacobi Equations

We define a ?-causal discretization of static convex Hamilton-Jacobi Partial Differential Equations (HJ PDEs) such that the solution value at a grid node is dependent only on solution values that are smaller by at least ?. We develop a Monotone Acceptance Ordered Upwind Method (MAOUM) that first determines a consistent, ?-causal stencil for each grid node and then solves the discrete equation in a single-pass through the nodes. MAOUM is suited to solving HJ PDEs efficiently on highly-nonuniform grids, since the stencil size adjusts to the level of grid refinement. MAOUM is a Dijkstra-like algorithm that computes the solution in increasing value order by using a heap to sort proposed node values. If ?>0, MAOUM can be converted to a Dial-like algorithm that sorts and accepts values using buckets of width ?. We present three hierarchical criteria for ?-causality of a node value update from a simplex of nodes in the stencil. The asymptotic complexity of MAOUM is found to be $\mathcal {O}((\hat{\Psi}\rho )^{d} N \log N)$ , where d is the dimension, $\hat{\Psi}$ is a measure of anisotropy in the equation, and ? is a measure of the degree of nonuniformity in the grid. This complexity is a constant factor $(\hat{\Psi}\rho)^{d}$ greater than that of the Dijkstra-like Fast Marching Method, but MAOUM solves a much more general class of static HJ PDEs. Although ? factors into the asymptotic complexity, experiments demonstrate that grid nonuniformity does not have a large effect on the computational cost of MAOUM in practice. Our experiments indicate that, due to the stencil initialization overhead, MAOUM performs similarly or slightly worse than the comparable Ordered Upwind Method presented in (Sethian and Vladimirsky, SIAM J. Numer. Anal. 41:323, 2003) for two examples on uniform meshes, but considerably better for an example with rectangular speed profile and significant grid refinement around nonsmooth parts of the solution. We test MAOUM on a diverse set of examples, including seismic wavefront propagation and robotic navigation with wind and obstacles.

Towards a Multi-Level Cache Performance Model for 3D Stencil Computation

It is crucial to optimize stencil computations since they are the core (and most computational demanding segment) of many Scientific Computing applications, therefore reducing overall execution time. This is not a simple task, actually it is lengthy and tedious. It is lengthy because the large number of stencil optimizations combinations to test, which might consume days of computing time, and the process is tedious due to the slightly different versions of code to implement. Alternatively, models that predict performance can be built without any actual stencil execution, thus reducing the cumbersome optimization task. Previous works have proposed cache misses and execution time models for specific stencil optimizations. Furthermore, most of them have been designed for 2D datasets or stencil sizes that only suit low order numerical schemes. We propose a flexible and accurate model for a wide range of stencil sizes up to high order schemes, that captures the behavior of 3D stencil computations using platform parameters. The model has been tested in a group of representative hardware architectures, using realistic dataset sizes. Our model predicts successfully stencil execution times and cache misses. However, predictions accuracy depends on the platform, for instance on x86 architectures prediction errors ranges between 1-20%. Therefore, the model is reliable and can help to speed up the stencil computation optimization process. To that end, other stencil optimization techniques can be added to this model, thus essentially providing a framework which covers most of the state-of-the-art.

A solution-adaptive lattice Boltzmann method for two-dimensional incompressible viscous flows

A stencil adaptive lattice Boltzmann method (LBM) is developed in this paper. It incorporates the stencil adaptive algorithm developed by Ding and Shu [26] for the solution of Navier–Stokes (N–S) equations into the LBM calculation. Based on the uniform mesh, the stencil adaptive algorithm refines the mesh by two types of 5-points symmetric stencils, which are used in an alternating sequence for increased refinement levels. The two types of symmetric stencils can be easily combined to form a 9-points symmetric structure. Using the one-dimensional second-order interpolation recently developed by Wu and Shu [27] along the straight line and the D2Q9 model, the adaptive LBM calculation can be effectively carried out. Note that the interpolation coefficients are only related to the lattice velocity and stencil size. Hence, the simplicity of LBM is not broken down and the accuracy is maintained. Due to the use of adaptive technique, much less mesh points are required in the simulation as compared to the standard LBM. As a consequence, the computational efficiency is greatly enhanced. The numerical simulation of two dimensional lid-driven cavity flows is carried out. Accurate results and improved efficiency are reached. In addition, the steady and unsteady flows over a circular cylinder are simulated to demonstrate the capability of proposed method for handling problems with curved boundaries. The obtained results compare well with data in the literature.

Trefftz difference schemes on irregular stencils

The recently developed Flexible Local Approximation MEthod (FLAME) produces accurate difference schemes by replacing the usual Taylor expansion with Trefftz functions -- local solutions of the underlying differential equation. This paper advances and casts in a general form a significant modification of FLAME proposed recently by Pinheiro & Webb: a least-squares fit instead of the exact match of the approximate solution at the stencil nodes. As a consequence of that, FLAME schemes can now be generated on irregular stencils with the number of nodes substantially greater than the number of approximating functions. The accuracy of the method is preserved but its robustness is improved. For demonstration, the paper presents a number of numerical examples in 2D and 3D: electrostatic (magnetostatic) particle interactions, scattering of electromagnetic (acoustic) waves, and wave propagation in a photonic crystal. The examples explore the role of the grid and stencil size, of the number of approximating functions, and of the irregularity of the stencils.

Treatment of Grid-Conforming Dielectric Interfaces in FDTD Methods

A new solution to the problem of implementing extended-stencil finite-difference operators across material interfaces is presented in this paper for time-domain electromagnetic simulations. In essence, modified spatial approximations are applied within a transition layer that encloses the discontinuity, while preserving the stencil size of the scheme employed in the homogeneous areas. The necessary difference expressions are determined with the aid of a straightforward design procedure, which is based on the reliable modeling of exact solutions, following a one-dimensional formulation. The proposed approach, whose beneficial effects are exemplified via numerical simulations, is analyzed for the case of the (2,4) method. Yet, its concepts are global and easily adapted to any discretization algorithm.

A Comprehensive New Methodology for Formulating FDTD Schemes With Controlled Order of Accuracy and Dispersion

Numerical dispersion errors in the wave-equation-finite-difference-time-domain (WE-FDTD) method have been treated by higher order schemes, coefficient modification schemes, dispersion relation preserving and non-standard schemes. In this work, a unified methodology is formulated for the systematic generation of WE-FDTD schemes tailored to the spectrum of the excitation. The methodology enables the scheme designer to gradually trade order of accuracy (OoA) for lower dispersion errors in a controlled manner at the cost of sacrificing low frequency behavior, that is not deemed critical for this type of excitation. The methodology is shown to encompass both existing and new schemes. Stability analysis is carried out concurrently with the generation of each scheme. Using a stencil size of 3 and 5 temporal and spatial samples, respectively, long term errors of a scheme designed for a specific pulse are compared with the standard (4,4) scheme that has the same computational complexity, via simulation of a modulated pulse that propagates over a million time steps.

A high-order discontinuous Galerkin method for the unsteady incompressible Navier–Stokes equations

We present a high-order discontinuous Galerkin discretization of the unsteady incompressible Navier–Stokes equations in convection-dominated flows using triangular and tetrahedral meshes. The scheme is based on a semi-explicit temporal discretization with explicit treatment of the nonlinear term and implicit treatment of the Stokes operator. The nonlinear term is discretized in divergence form by using the local Lax–Friedrichs fluxes; thus, local conservativity is inherent. Spatial discretization of the Stokes operator has employed both equal-order ( P k − P k ) and mixed-order ( P k − P k−1 ) velocity and pressure approximations. A second-order approximate algebraic splitting is used to decouple the velocity and pressure calculations leading to an algebraic Helmholtz equation for each component of the velocity and a consistent Poisson equation for the pressure. The consistent Poisson operator is replaced by an equivalent (in stability and convergence) operator, namely that arising from the interior penalty discretization of the standard Poisson operator with appropriate boundary conditions. This yields a simpler and more efficient method, characterized by a compact stencil size. We show the temporal and spatial behavior of the method by solving some popular benchmarking tests. For an unsteady Stokes problem, second-order temporal convergence is obtained, while for the Taylor vortex test problem on both semi-structured and fully unstructured triangular meshes, spectral convergence with respect to the polynomial degree k is obtained. By studying the Orr–Sommerfeld stability problem, we demonstrate that the P k − P k method yields a stable solution, while the P k − P k−1 formulation leads to unphysical instability. The good performance of the method is further shown by simulating vortex shedding in flow past a square cylinder. We conclude that the proposed discontinuous Galerkin method with the P k − P k formulation is an efficient scheme for accurate and stable solution of the unsteady Navier–Stokes equations in convection-dominated flows.

Fourth order partial differential equations on general geometries

We extend a recently introduced method for numerically solving partial differential equations on implicit surfaces [M. Bertalmío, L.T. Cheng, S. Osher, G. Sapiro. Variational problems and partial differential equations on implicit surfaces, J. Comput. Phys. 174 (2) (2001) 759–780] to fourth order PDEs including the Cahn–Hilliard equation and a lubrication model for curved surfaces. By representing a surface in RN as the level set of a smooth function, ϕ, we compute the PDE using only finite differences on a standard Cartesian mesh in RN. The higher order equations introduce a number of challenges that are of less concern when applying this method to first and second order PDEs. Many of these problems, such as time-stepping restrictions and large stencil sizes, are shared by standard fourth order equations in Euclidean domains, but others are caused by the extreme degeneracy of the PDEs that result from this method and the general geometry. We approach these difficulties by applying convexity splitting methods, ADI schemes, and iterative solvers. We discuss in detail the differences between computing these fourth order equations and computing the first and second order PDEs considered in earlier work. We explicitly derive schemes for the linear fourth order diffusion, the Cahn–Hilliard equation for phase transition in a binary alloy, and surface tension driven flows on complex geometries. Numerical examples validating our methods are presented for these flows for data on general surfaces.

A convergent monotone difference scheme for motion of level sets by mean curvature

An explicit convergent finite difference scheme for motion of level sets by mean curvature is presented. The scheme is defined on a cartesian grid, using neighbors arranged approximately in a circle. The accuracy of the scheme, which depends on the radius of the circle, dx, and on the angular resolution, dθ, is formally O(dx2+dθ). The scheme is explicit and nonlinear: the update involves computing the median of the values at the neighboring grid points. Numerical results suggest that despite the low accuracy, acceptable results are achieved for small stencil sizes. A numerical example is presented which shows that the centered difference scheme is non-convergent.

Refinable bivariate quartic $C^2$-splines for multi-level data representation and surface display

In this paper, a second-order Hermite basis of the space of C 2 -quartic splines on the six-directional mesh is constructed and the refinable mask of the basis functions is derived. In addition, the extra parameters of this basis are modified to extend the Hermite interpolating property at the integer lattices by including Lagrange interpolation at the half integers as well. We also formulate a compactly supported super function in terms of the basis functions to facilitate the construction of quasi-interpolants to achieve the highest (i.e., fifth) order of approximation in an efficient way. Due to the small (minimum) support of the basis functions, the refinable mask immediately yields (up to) four-point matrix-valued coefficient stencils of a vector subdivision scheme for efficient display of C 2 -quartic spline surfaces. Finally, this vector subdivision approach is further modified to reduce the size of the coefficient stencils to two-point templates while maintaining the second-order Hermite interpolating property.

Comparison of the Dispersion Properties of Higher Order FDTD Schemes and Equivalent-Sized MRTD Schemes

The dispersion errors of higher order finite-difference time-domain (HO-FDTD) algorithms are compared to those of multiresolution time-domain (MRTD) algorithms that have equivalent spatial stencil sizes. Both scaling-function-based MRTD (S-MRTD) and wavelet-function-based MRTD (W-MRTD) schemes are considered. In particular, the MRTD schemes considered include the Coifman scaling functions and the Cohen-Daubechies-Feauveau (CDF) biorthogonal scaling and wavelet functions. In general, the HO-FDTD schemes are more accurate than their MRTD counterparts.

Finite-volume compact schemes on staggered grids

Compact finite-difference schemes have been recently used in several Direct Numerical Simulations of turbulent flows, since they can achieve high-order accuracy and high resolution without exceedingly increasing the size of the computational stencil. The development of compact finite-volume schemes is more involved, due to the appearance of surface and volume integrals. While Pereira et al. [J. Comput. Phys. 167 (2001)] and Smirnov et al. [AIAA Paper, 2546, 2001] focused on collocated grids, in this paper we use the staggered grid arrangement. Compact schemes can be tuned to achieve very high resolution for a given formal order of accuracy. We develop and test high-resolution schemes by following a procedure proposed by Lele [J. Comput. Phys. 103 (1992)] which, to the best of our knowledge, has not yet been applied to compact finite-volume methods on staggered grids. Results from several one- and two-dimensional simulations for the scalar transport and Navier–Stokes equations are presented, showing that the proposed method is capable to accurately reproduce complex steady and unsteady flows.

A New Compact Spectral Scheme for Turbulence Simulations

We propose a new kind of compact difference scheme for the computation of the first and second derivatives in the simulation of high-Reynolds number turbulent flows. The scheme combines and truncates the pseudospectral representation of derivative for convergence acceleration. Comparison of the wave resolution property with available optimized compact schemes minimizing the prescribed wave resolution error reveals our scheme's superiority for the same size of stencils without introducing optimization parameters. An accompanying boundary scheme is also proposed with the stability analysis. The proposed scheme is tested for the evaluation of derivatives of a function that decays very slowly in the wavenumber space, and for the simulation of three-dimensional isotropic turbulence.

Prefactored Small-Stencil Compact Schemes

A new type of prefactored compact schemes is presented that obtains high-order accuracy while using a very small stencil size. These schemes require fewer boundary stencils and offer simpler boundary condition implementation than the existing compact schemes. Results are shown for linear and nonlinear 1-D benchmark problems and a 2-D linear benchmark problem. Boundary stencils for both exterior and interior multiblock boundaries are given, and boundary condition specification is described and demonstrated.

On numerical techniques in CFD

Numerical techniques play an important role in CFD. Some of them are reviewed in this paper. The necessity of using high order difference scheme is demonstrated for the study of high Reynolds number viscous flow. Physical guide lines are provided for the construction of these high order schemes. To avoid unduly ad hoc tremtment in the boundary region the use of compact scheme is recommended because it has a small stencil size compared with the traditional finite difference scheme. Besides preliminary Fourier analysis shows the compact scheme can also yield better space resolution which makes it more suitable to study flow with multiscales e.g. turbulence. Other approaches such as perturbation method and finite spectral method are also emphasized. Typical numerical simulations were carried out. The first deals with Euler equations to show its capabilities to capture flow discontinuity. The second deals with Navier-Stokes equations studying the evolution of a mixing layer, the pertinent structures at different times are shown. Asymmetric break down occurs and also the appearance of small vortices.

Ion projection lithography for resistless patterning of thin magnetic films

The ion projector at the Fraunhofer Institute in Berlin (manufactured by Ion Microfabrication Systems, IMS, Vienna) has been used to pattern thin magnetic films. Disordering of chemically ordered FePt films by patterned ion irradiation (dose: 10 16 He + ions/cm 2; 75 keV) without a resist mask on the sample is demonstrated. Magnetically altered areas have been investigated by magnetic force microscopy. Magnetic islands of 340 nm width could be detected, limited by the size of the open stencil mask holes. The projector is capable of 50 nm resolution which would result in a storage density of 64 Gb/in 2. In a second experiment the ion projector has been applied for structured ion milling with Xe + ions on thin metallic films also without a resist mask. In a polycrystalline Au film 130 nm lines and spaces with a depth of 8 nm could be produced by a dose of 10 15Xe +ions at 75 keV energy.

A Wavenumber Based Extrapolation and Interpolation Method for Use in Conjunction with High-Order Finite Difference Schemes

The errors incurred in using extrapolation and interpolation in large scale computations are analyzed and quantified in the wavenumber space. If a large extrapolation stencil is used, the errors in the low wavenumbers can be significantly reduced. However, the errors in the high wavenumbers are, at the same time, greatly increased. The opposite is true if the stencil size is reduced. Based on the wavenumber analysis, an optimized extrapolation and interpolation method is proposed. The optimization is carried out over a selected band of wavenumbers. It is known that extrapolation often leads to numerical instability. The instability is the result of large error amplification in the high wavenumber range. To reduce the tendency to trigger numerical instability, it is proposed that an extra constraint be imposed on the optimized extrapolation method. The added constraint aims to reduce error amplification over the high wavenumbers. Numerical examples are provided to illustrate that accurate and stable numerical results can be obtained in large scale simulation using a high-order finite difference scheme and the proposed optimized extrapolation method. When the same problems are recomputed using the familiar high-order polynomials extrapolation method in the Lagrange form, in one case the numerical results are plagued by large errors and ultimately instability. In another problem, it is found that the use of the Lagrange polynomials extrapolation method would lead immediately to numerical instability.

Stability and dispersion analysis of Battle-Lemarie-based MRTD schemes

The stability and dispersion performance of the recently developed Battle-Lemarie multiresolution time-domain schemes is investigated for different stencil sizes. The contribution of wavelets is enhanced and analytical expressions for the maximum allowable time step are derived. It is observed that larger stencils decrease the numerical phase error, making it significantly lower than finite-difference time domain for low and medium discretizations. The addition of wavelets further improves the dispersion performance for discretizations close to the Nyquist limit, though it decreases the value of the maximum time step, guaranteeing the stability of the scheme.

An adaptive multigrid technique for evaluating long-range forces in biomolecular simulations

We propose an adaptive multigrid method for evaluating the nonlocal, long-range interactions in biomolecular models. The method provides O( h 3) accuracy for the forces, where h is the mesh size of the finest cubic grid, even though the stencil size is minimal. Two new algorithmic components are developed: an adaptive technique for spreading a three-dimensional charge distribution onto a space mesh, and a procedure for performing coarse-to-fine grid interpolation of the potential, force components, and their derivatives (grid-to-particle interpolation at the finest level) in a consistent fashion. The adaptive mesh used in the spreading procedure contains two types of nodes: regular, positioned at nodes cubic grid, and flexible, one per cubic cell. Thus, the charges located in each grid cell are spread to these nine nodes. The method conserves charge, dipole, and quadrupole moments of cells. The spreading and interpolation techniques are economical in computation time and storage due to the overlap design in the stencils. Implementation is also straightforward. Our initial tests demonstrate the favorable performance for low accuracy with respect to the fast multipole technique.