A parallel geometric multigrid method for finite elements on octree meshes applied to elastic image registration

The first component of this thesis is a parallel algorithm for constructing octree meshes for finite element computations. Prior to octree meshing, the linear octree data structure must be constructed and a constraint known as “2:1 balancing” must be enforced; parallel algorithms for these two subproblems are also presented. The second component of this thesis is a parallel geometric multigrid algorithm for solving elliptic partial differential equations (PDEs) using these octree meshes. The last component of this thesis is a parallel multiscale Gauss Newton optimization algorithm for solving the elastic image registration problem. The registration problem is discretized using finite elements on octree meshes and the parallel geometric multigrid algorithm is used as a preconditioner in the Conjugate Gradient (CG) algorithm to solve the linear system of equations formed in each Gauss Newton iteration. The parallel octree meshing and multigrid algorithms have several physical and computer science applications such as in solid/fluid mechanics, heat/mass transfer, electromagnetism, image processing and unstructured mesh generation. Potential applications for the image registration algorithm include automatic identification of abnormalities in medical images, motion reconstruction from temporal sequences of images and planning of surgeries. Several ideas were used to reduce the overhead for constructing the octree meshes. These include (a) a way to lower communication costs by reducing the number of synchronizations and reducing the communication message size, (b) a way to reduce the number of searches required to build element-to-vertex mappings, and (c) a compression scheme to reduce the memory footprint of the entire data structure. To our knowledge, the multigrid algorithm presented in this work is the only matrix-free multiplicative geometric multigrid implementation for solving finite element equations on octree meshes using thousands of processors. The overall scheme is second-order accurate, for sufficiently smooth right-hand sides and material properties; and its complexity, for nearly uniform trees, is OpNnplog Nnpp+O (np log np). Here, N is the number of octants and np is the number of processors. The proposed registration algorithm is also unique; it is a combination of many different ideas: adaptivity, parallelism, fast optimization algorithms, and fast linear solvers. All the algorithms were implemented in C++ using the Message Passing Interface (MPI) standard and were built on top of the PETSc library from Argonne National Laboratory. The multigrid implementation has been released as an open source software: Dendro. Several numerical experiments were performed to test the performance of the algorithms. These experiments were performed on a variety of NSF TeraGrid platforms: on the Cray XT3 MPP system “Bigben” at the Pittsburgh Supercomputing Center (PSC), the Intel 64 Linux Cluster “Abe” at the National Center for Supercomputing Applications (NCSA), and the Sun Constellation Linux Cluster “Ranger” at the Texas Advanced Computing Center (TACC). Our largest run was a highly-nonuniform, 8-billion-unknown, elasticity calculation on 32,000 processors.