State of the art Stokes Solver using FMM

Stokes equation is one of the most important equations derived from Navier-Stokes. Numerical solutions and descritazation of the Stokes equation is challenging. For instance, one cannot use arbitrary descretization spaces for velocity and pressure. It is an elliptic but indefinite problem, which further complicates the construction of fast linear algebraic solvers and preconditioners, especially for problems with highly variable coefficients or high-order discretizations. We are developing an integral formulation that circumvents most of the difficulties with the Stokes equation. This is done by using pvfmm, a novel adaptive volume fast multipole method. Compared to finite element methods, our formulation decouples the velocity and pressure, and generates fields that are by construction divergence free. Our benchmarking tests for Stokes problem have shown a four fold speed up compared to deal.II, a very well written and leading FEM software.

To illustrate the capabilities of our method, we simulated Stokes flow through a porous medium. Figure 1 shows the porous medium geometry used. The streamlines of the flow are shown in Figure 2, and the mesh structure and the parallel partitioning are shown in Figure 3.

Update: Our recent paper submitted to SC14 has been accepted and nominated as a finalist for best student paper. You can download a draft of the paper here.

Close Project

Brain Tumor Inverse Problem

Gliomas are tumors that arise from Glial cells in the brain. They accounted for 29% of all brain and central nervous system (CNS) tumors, and 80% of all malignant tumors out of about 60,000 cases diagnosed each year in the US. Gliomas are classified into four grades by the World Health Organization (WHO). Grade I represents the most benign and curable Glioma. However, grade II-IV are currently not curable, in part because of their aggressive infiltration into normal (healthy) tissue. Despite advances in surgery, chemo/radio therapy, the median survival rate of high grade Gliomas has remained about one year in the past 30 years.

One of the key parameters in increasing the survival rate of patients diagnosed with malignant manifestations of Glioma is how well the tumor invasion boundaries are detectable. Unfortunately, with the current imaging technologies only the bulk of the tumor abnormalities, can be detected, and the infiltrated tumor cells get masked. These are the cells that cause tumor recurrence after treatment. I am trying to address the problem of determining/approximating the full extent of the tumor infiltration using the data available from different imaging modalities, as well as the tumor growth dynamics.

• Our paper has recently been accepted in the Journal of Mathematical Biology. Please click here to download the paper.

Figure 2: Simulation of tumor growth at three different time frames. The rows show tumor distribution at t=0, 1, and 2 (or in dimensional form 0, 14, and 28 months, respectively).

Close Project

AccFFT: A New Massively Parallel FFT Library for CPU and GPU Architectures

Accelerated FFT (AccFFT) is a new massively parallel FFT library for computing distributed Fast Fourier Transforms. The library has been designed with the goal of achieving maximum performance, without making the user interface complicated. AccFFT supports parallel FFTs distributed with slab or pencil decomposition for both CPU and GPU architectures. The library's scalability has been tested upto 131K cpu cores, and upto 4K GPUs. It achieves a speedup of upto 4 times compared to other implementations.

Close Project

Parallel Image Registration

Figure 1: (left) Visualization of image registration problem. The goal is to find y, which is a mapping from template image ρT to the reference image ρR. (right): Parallel decomposition of the domain is shown along with a single step trajectory of our Semi-Lagrangian solver. The three trajectories show how the grid deformed at a prior time step. Except for the black trajectory, a sparse all-to-all communication is needed to compute the results, since the deformations can cross boundaries of the processors (e.g. the red trajectory).

Image registration is a process in which a mapping from a reference image to a target image is sought. For instance, one may wish to align a set of images as a preprocessing step for a machine learning task, or to align MRI images of a patient taken at different times. In certain cases image registration can also be used to deduce missing properties in an image such as in tumorous parts, by registering the patient’s brain to a healthy atlas. We formulate the problem as a PDE constrained optimization. We treat the problem as finding a variable velocity that can advect the template image to the reference image. Only diffeomorphic transformation are allowed, which requires positive determinant of the deformation map. Solving the image registration problem in 3D is a very time consuming task with a prohibitive time to solution. We use second order reduced Hessian method to solve the optimization problem. Using the Hessian leads to second order convergence in the vicinity of the MAP point, however it requires expensive forward and adjoint solves in time. As discussed above, we use both mathematics and computer science to address this. For the former we use novel Hessian preconditioners to increase the convergence rate of each Hessian solve. For the latter, we have developed a novel parallel Semi Lagrangian scheme that is unconditionally stable with excellent parallel scalability.
Semi Lagrangian is a variant of particle methods in which the properties are found by solving for the trajectory backwards in time, as shown schematically in Fig. 1. In the case of distributed-memory parallelism, each processor only has access to disjoint parts of the domain. Therefore, it is possible that the trajectory would fall outside its locally owned space. To solve the advection problem one needs to perform parallel interpolation to get the desired values (for instance the intensity of the image at current time or its deformation). The code that we have developed is based on AccFFT along with a parallel high-order interpolation. The latter uses sparse alltoall collective to compute the interpolated values that map outside each processor’s domain. Using this approach, we have solved a registration problem of size 256³ in less than five seconds on 32 nodes of Maverick system at TACC, which is quite remarkable. The result of this work was presented in this year’s SuperComputing conference [pdf].

Close Project

A Comparative Study of Massively Parallel Poisson Solvers

The need for large scale parallel solvers for elliptic partial differential equations (PDES) pervades across a spectrum of problems with resolution requirements that cannot be accommodated on current systems. Poisson solvers must scale to trillions of unknowns. Example of methods that scale well are the FFT (based on spectral discretizations), the Fast Multipole Method, and multigrid methods (for stencil-based discretizations). We have benchmarked these methods and compared their parallel efficiency as well as the corresponding cost per unknowns for different test cases. FFT is tested with AccFFT, FMM with pvfmm, AMG with ML package, and GMG with an in house code. The results of the comparison can be found in this paper (submitted to SISC).

Close Project