Multigrid methods and optimization
In many situations in scientific computing, it is often necessary to solve sparse linear systems that are represented by a matrix equation Ax = b, where x is an unknown vector of length N, b is a known vector of length N, and A is an N × N matrix where most entries are zero. An example where this would arise would be in electrostatics, where the vector b would represent a distribution of charges and the vector x would represent the electric potential.
There are many methods of solving systems such as these. Simple methods, such as the Jacobi or Gauss–Seidel methods, start with an initial guess for x and then slowly improve it via iteration. However, while these methods are able to rapidly remove high frequency errors in the solution, they take a long time to remove low frequency errors, requiring O(N) computation time per grid point. I have been particularly interested in the use of the multigrid algorithm, which often provides highly efficient solution of these systems and is particularly well-suited to many physical problems. By employing a hierarchy of progressively coarser grids, it is able to damp both high frequency and low frequency errors effectively, and can result in O(1) computation time per grid point. Below, I describe several research projects that I have carried out that rely on the multigrid algorithm.
A multigrid algorithm on adaptive quadtree and octree grids
In many physical simulation problems, it is advantageous to use adaptive grids, where the resolution is increased in some particular areas of interest—this can often result in dramatic savings in the amount of computational effort required to solve a problem. One method of creating an adaptive grid is to start with a square, and the repeatedly subdivide the square into groups of four squares. An example grid is shown here, where the resolution of the grid has been chosen to be higher near a boundary of interest, shown in red. Grids such as these can be efficiently represented on a computer via the use of a quadtree data structure. The same principle can be applied in three dimensions by starting with a cube that can be successively subdivided into groups of eight cubes, and using an octree data structure.
In collaboration with Frédéric Gibou and Maxime Theillard at UC Santa Barbara, I worked on an implementation of the multigrid algorithm on grids such as these. The quadtree and octree data structures make it very easy to define a hierarchy of grids for use with the multigrid algorithm, and in a recent paper, we were able to demonstrate rapid convergence for a variety of test problems in two and three dimensions.
- M. Theillard, C. H. Rycroft, and F. Gibou, A multigrid method on non-graded adaptive octree/quadtree cartesian grids, J. Sci. Comput. 55, 1–15 (2013). [Link]
Computation of three-dimensional time-periodic water waves
Water waves in the ocean have a surprisingly complex structure, arising from nonlinearities in the underlying equations of motion, and they have attracted much interest from mathematicians and scientists. One area of study has been to search for perfectly time-periodic waves, and in 1880, George Gabriel Stokes proposed that the maximum time-periodic traveling wave would have a crest with an internal angle of 120°, a result that was later proved analytically and investigated numerically. There have been many other numerical studies of time-periodic wave motion, although most of these have focused on the two-dimensional case, and there are no examples of large-amplitude standing waves being computed in three dimensions, mainly due to the high computational cost of carrying out such a search.
I have worked with Jon Wilkening to address this problem, by creating some new numerical algorithms. In our simulations, we make track the evolution of the wave height and velocity. To time-evolve the wave, we make use of spectral methods, but also need to compute a velocity potential in the bulk, which is done using a fourth-order curved finite-element method. This finite-element problem is very large, and is represented by a sparse linear system with up to a billion non-zero entries. To solve this, I wrote a parallel geometric multigrid algorithm in C++ using the OpenMPI library. Levenberg–Marquardt optimization methods were then employed to search for several families of three-dimensional standing waves, highlighting a number of physical effects such as resonances, and notable differences with the two-dimensional counterparts. We are currently extending this methodology to examine more general questions of wave motion, and are making use of similar techniques to find time-periodic solutions of the compressible Euler equations.