Adaptive Mesh Refinement for Characteristic Grids

I consider techniques for Berger-Oliger adaptive mesh refinement (AMR) when numerically solving partial differential equations with wave-like solutions, using characteristic (double-null) grids. Such AMR algorithms are naturally recursive, and the best-known past Berger-Oliger characteristic AMR algorithm, that of Pretorius & Lehner (J. Comp. Phys. 198 (2004), 10), recurses on individual "diamond" characteristic grid cells. This leads to the use of fine-grained memory management, with individual grid cells kept in 2-dimensional linked lists at each refinement level. This complicates the implementation and adds overhead in both space and time. Here I describe a Berger-Oliger characteristic AMR algorithm which instead recurses on null \emph{slices}. This algorithm is very similar to the usual Cauchy Berger-Oliger algorithm, and uses relatively coarse-grained memory management, allowing entire null slices to be stored in contiguous arrays in memory. The algorithm is very efficient in both space and time. I describe discretizations yielding both 2nd and 4th order global accuracy. My code implementing the algorithm described here is included in the electronic supplementary materials accompanying this paper, and is freely available to other researchers under the terms of the GNU general public license.Comment: 37 pages, 15 figures (40 eps figure files, 8 of them color; all are viewable ok in black-and-white), 1 mpeg movie, uses Springer-Verlag svjour3 document class, includes C++ source code. Changes from v1: revised in response to referee comments: many references added, new figure added to better explain the algorithm, other small changes, C++ code updated to latest versio

On the error propagation of semi-Lagrange and Fourier methods for advection problems

In this paper we study the error propagation of numerical schemes for the advection equation in the case where high precision is desired. The numerical methods considered are based on the fast Fourier transform, polynomial interpolation (semi-Lagrangian methods using a Lagrange or spline interpolation), and a discontinuous Galerkin semi-Lagrangian approach (which is conservative and has to store more than a single value per cell). We demonstrate, by carrying out numerical experiments, that the worst case error estimates given in the literature provide a good explanation for the error propagation of the interpolation-based semi-Lagrangian methods. For the discontinuous Galerkin semi-Lagrangian method, however, we find that the characteristic property of semi-Lagrangian error estimates (namely the fact that the error increases proportionally to the number of time steps) is not observed. We provide an explanation for this behavior and conduct numerical simulations that corroborate the different qualitative features of the error in the two respective types of semi-Lagrangian methods. The method based on the fast Fourier transform is exact but, due to round-off errors, susceptible to a linear increase of the error in the number of time steps. We show how to modify the Cooley--Tukey algorithm in order to obtain an error growth that is proportional to the square root of the number of time steps. Finally, we show, for a simple model, that our conclusions hold true if the advection solver is used as part of a splitting scheme.Comment: submitted to Computers & Mathematics with Application

Numerical Methods for the Fractional Laplacian: a Finite Difference-quadrature Approach

The fractional Laplacian $(-\Delta)^{\alpha/2}$ is a non-local operator which depends on the parameter $\alpha$ and recovers the usual Laplacian as $\alpha \to 2$. A numerical method for the fractional Laplacian is proposed, based on the singular integral representation for the operator. The method combines finite difference with numerical quadrature, to obtain a discrete convolution operator with positive weights. The accuracy of the method is shown to be $O(h^{3-\alpha})$. Convergence of the method is proven. The treatment of far field boundary conditions using an asymptotic approximation to the integral is used to obtain an accurate method. Numerical experiments on known exact solutions validate the predicted convergence rates. Computational examples include exponentially and algebraically decaying solution with varying regularity. The generalization to nonlinear equations involving the operator is discussed: the obstacle problem for the fractional Laplacian is computed.Comment: 29 pages, 9 figure

How to mesh up Ewald sums (I): A theoretical and numerical comparison of various particle mesh routines

Standard Ewald sums, which calculate e.g. the electrostatic energy or the force in periodically closed systems of charged particles, can be efficiently speeded up by the use of the Fast Fourier Transformation (FFT). In this article we investigate three algorithms for the FFT-accelerated Ewald sum, which attracted a widespread attention, namely, the so-called particle-particle-particle-mesh (P3M), particle mesh Ewald (PME) and smooth PME method. We present a unified view of the underlying techniques and the various ingredients which comprise those routines. Additionally, we offer detailed accuracy measurements, which shed some light on the influence of several tuning parameters and also show that the existing methods -- although similar in spirit -- exhibit remarkable differences in accuracy. We propose combinations of the individual components, mostly relying on the P3M approach, which we regard as most flexible.Comment: 18 pages, 8 figures included, revtex styl