Research (This description is somewhat dated...)


My primary research interest is the development of fast, stable and accurate numerical algorithms for the approximation of partial differential equations arising in engineering and natural sciences. A large part of my research has been devoted to the development of artificial boundary conditions necessary for solution of time dependent partial differential equations on unbounded domains. I have also worked on several inter-disciplinary projects developing both numerical methods and physical models. Recently I have become interested in nucleation phenomena associated with phase transitions in complex fluids and biological systems that we study by combining numerical simulations of self-consistent field theory (SCFT) with minimum energy path (MEP) methods.

Arbitrary order Hermite methods

This project concerns the development of a massively parallel arbitrary-order subsonic compressible Navier-Stokes solver based on Hermite discretizations. The solver is currently used for direct simulations of turbulent jet noise at moderate Reynolds numbers. Such simulations are extremely challenging and must be performed on supercomputers. Hermite methods are ideal for exploiting present and future multicore supercomputer architectures as they posses the following key features:

  • High-resolution enabling a minimal number of degrees-of-freedom per wavelength;
  • Large computation-to-communication ratio to efficiently use multicore nodes and local memory;
  • Flexibility and modularity enabling architecture-specific optimization of core computational kernels for use in a wide variety of applications.

Being new and unexplored methods, there are many exciting aspects that need to be examined: how to incorporate shock capturing, adaptive mesh refinement, stability preserving time-stepping, non-reflecting boundary conditions, how to best handle geometry, how to establish non-linear stability. Finding satisfactory answers to all these questions will be a rewarding task spanning over the next few years. Some recent results, for example hybridization with nodal discontinuous Galerkin methods can be found on the the publications page.
Collaborators: Tom Hagstrom, SMU, Tim Colonius, Caltech.

Numerical simulations in self-consistent field theory and study of rare transitions in polymer systems

Human gene therapy shows broad potential as a primary means of treating disease. Polymer-based vectors rely on the ability to electrostatically bind and condense the genetic material into nanoparticles, which are termed ``polyplexes''. Designing an efficient gene-delivery vector requires a mechanistic understanding of each step during the gene delivery process. In particular, the nanoparticle must interact with, and ultimately rupture, the biomembrane before entering the nucleus. Accurate numerical SCF solutions can reveal the thermodynamic states of this system, but to fully understand the problem requires also knowing the energies of the transition states between these (meta)stable states. These high energy states cannot be obtained by the SCFT, but require additional methods, such as the string method, to obtain the so called minimum energy path. The string method relies on the evolution of states based on the gradient of the energy.
Collaborators: Christina Ting and Zhen-Gang Wang, Caltech.

Numerical methods for PDE and their application to engineering problems

Numerical methods for elastic waves

As a postdoctoral scholar at Lawrence Livermore National Laboratory I developed numerical methods and massively parallel software for seismic modeling. My main theoretical contribution to the project was the development of a bullet-proof explicit finite difference discretization for the elastic wave equation with free surfaces on curvilinear grids. Using summation by parts techniques we could prove that our fully discrete approximation preserved a discrete energy to machine precision for media with arbitrarily-varying physical properties. Stability for such media is critical in seismic applications where the material properties vary wildly. The offspring the open source massively parallel 3D code WPP, developed in collaboration with Anders Petersson and Björn Sjögreen, is routinely used on thousands of cpus by seismologists at LLNL.

I have also incorporated the stable discretization into an open-source 2D/3D parallel overset-grid-solver for solid mechanics computations on complex geometries. The software is written within the Overture framework. This work was performed in collaboration with Bill Henshaw, Jeff Banks and Don Schwendeman.

Compact Difference Embedded Boundary Methods

Finite difference discretizations with the boundary embedded in a Cartesian grid provide an alternative to body-fitted overset-grid methods, finite-element and finite-volume methods for handling complex geometries. Embedded boundary methods, being based on structured Cartesian grids, are highly efficient both in terms of operations and required memory per degree of freedom. They are also well suited for massively parallel computers as they are easy to load balance and scale to thousands of cores. Another advantage is that the geometry can be represented in a lower dimension and only local information such as the location and the normal of the boundary are required to enforce the boundary conditions, no costly (parallel) grid generation is needed.

We have developed a high-order-accurate embedded boundary method based on compact finite-difference approximations of spatial derivatives. We use two different approaches, related to the second-order-accurate methods described by Kreiss, Petersson and Yström. One notable difference is that we don't assign the ghost values by extrapolation to the outside of the boundary, but rather impose boundary conditions by assigning values via interpolation to "ghost points" just inside the boundary. Our new approach is more accurate and removes the small cell stiffness problem for the Dirichlet problem as the boundary is always well separated from the point interior to the ghost point. For the Neumann problem we mitigate the small cell problem by adding a fifth order accurate regularizing term to the interpolant before taking the normal derivative of the interpolant. An advantage with our method is that both Neumann and Dirichlet boundary conditions can be handled within the same framework.

Converging shock waves

Converging shock waves is an effective method to generate high temperatures and pressures for experimental, medical and engineering purposes. To maintain a predictable shape of the shock during the focusing event it is important that the shape is stable. With Eliasson and Henshaw we studied, numerically, the stability of converging polygonally shaped shocks by solving the Euler equations of gas dynamics using a high-order accurate Godunov method. The geometry was discretized with overlapping structured grids and adaptive mesh refinement was used to dynamically track the shocks and contact surfaces. Results from the simulations are shown below for an initially cylindrical shock diffracted by 3, 8 and 16 cylinders.

Computational modeling of natural convection for a Titan Montgolfiere balloon

In a joint project with the Jet Propulsion Laboratory we developed a computational model to predict the buoyant flow induced inside and around balloons with idealized heat sources. Our model will be used to verify existing JPL system level models and potentially as a tool to understand the flow physics of different balloon designs; for example single- vs. double-walled balloons. The full governing equations describing conservation of mass, momentum, and energy are solved internal and external to the balloon membrane. We consider a wide range of heat sources that include both laminar and turbulent flows.

Truncation of unbounded domains using absorbing layers

Many problems in physics, biology and other natural sciences are described by equations posed on unbounded domains. To compute a numerical solution to such problems it is necessary, due to finite computational resources, to truncate the unbounded domain. A popular approach is to surround the truncated domain with an absorbing layer where the solution decays rapidly so that reflections from the outermost boundary are negligible.

Perfectly Matched Layers: General formulation, well-posedness and stability

The, by far, most popular absorbing layer is Berenger's Perfectly Matched Layer (PML) for electro-magnetics. The perfect matching guarantees that the waves propagate without reflection through the interface between the PML and the computational domain. A significant part of my research has focused on the generalization and application of the PML method. The overarching goal of this research has been to answer the fundamental question: For which problems can a well-posed, stable and efficient PML be constructed?

To study this question a very general PML model, including arbitrary many free parameters was used. We found that by including a certain parameter in the direction tangential to the layer it could be proved that the PML equations are strongly well-posed as long as the original hyperbolic system is also strongly well-posed. This result concluded a long standing discussion of weak and strong well-posedness of general PML models initiated a decade ago by Abarbanel and Gottlieb.

We also considered the stability of our PML model. Obtaining stability is a delicate process and, in general, each application must be examined separately. To simplify such investigations we developed a method for automatically symmetrizing Petrowskii well-posed Cauchy problems. The method is rooted in the Sturm sequence technique for determining the location of the roots of a complex polynomial and can also be used to decide stability (e.g. of a PML) by verifying a small number of algebraic inequalities. If the inequalities hold, they define coefficients in a local energy density from which an energy estimate can be automatically constructed. The method is easy to apply and we have used it to prove stability for PML applied to: Maxwell's equations, the linearized Euler equations and arbitrary 2 by 2 symmetric hyperbolic systems in (2+1) dimensions.

The most recent generalization of our PML model is to hyperbolic-parabolic systems. This is important as it allows for the treatment of the compressible Navier-Stokes equations. We have also studied PML for non-linear dispersive wave equations, and the elastic wave equation.

The high-order super-grid-scale method

An alternative to PML is the super-grid-scale method introduced by Colonius and Ran. In the super-grid-scale method the domain truncation problem is treated as a modeling problem rather than a mathematical limit process. Its main advantage is that the super-grid-scale method is straightforward to use for applications where no exact mathematical boundary condition is available. This is the case in most applications where the physical properties, models and governing equations are nonlinear, change or are under development. As an extension of the original super-grid we design a well-posed and stable high-order-accurate super-grid-scale method for linear hyperbolic systems. We show that a subset of the model is exact for the continuous problem and that the polluting super-scales arise from the discretization itself. To mitigate this we introduced a modeling term, consisting of a high-order artificial viscosity, tailored to selectively damp unresolved waves. With this simple yet powerful approach we achieve PML like accuracy for problems where PML is unstable. In particular we demonstrate that the suggested method converge at the order of the scheme used in the interior