High Order Spatial Generalization of 2D and 3D Isotropic Discrete Gradient Operators with Fast Evaluation on GPUs

Based on the concept of isotropic centered finite differences, this work generalizes the spatial order of accuracy of the 2D and 3D isotropic discrete gradient operators to a higher order. A suitable methodology is used to obtain a set of equations from which it is possible to deduce stencil weights to achieve numerical approximations of both high order spatial and high order isotropic gradients. We consider that the suggested discretization will be useful for enhancing the quality of the results in various scientific fields. The spatial order (S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S$$\end{document}) controls the spatial order of accuracy of the gradient norm and direction, while the isotropic order (I\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I$$\end{document}) controls, in some situations, the spatial order of accuracy of the gradient direction. A useful list of the stencil weights needed to construct different high order spatial and isotropic gradients is given. Numerical tests show that the numerical spatial orders of accuracy of the gradient approximation are the same as those predicted theoretically. Also, to illustrate the benefit of the new discretizations, some simulations with a multiphase lattice Boltzmann model are presented. Then, a series of benchmarks comparing various efficient convolution algorithms used to compute function or image gradients is presented. Different platforms implemented on CPUs and GPUs are studied, namely: plain MATLAB; the Jacket plugin for MATLAB; and CUDA. The results show situations in which substantial computational speedup can be obtained with CUDA and the Jacket plugin for MATLAB versus MATLAB on a CPU. Examples of 2D and 3D gradient computations using convolution products performed with our code are available for download as electronic supplementary material.