I am a Senior Research Scientist at Pixar Research.
I research topics in physics-based simulation, including fire, water, and humans.
My work has appeared in over two dozen movies, and received a 2012 SciTech Oscar. Many of the results can be seen on YouTube,
and extensive source code is available on this page.

Previously, I was an Associate Professor with tenure at the University of California at Santa Barbara,
where I was affiliated with the Media Arts and Technology Program and the Department of Computer Science.
While there, I received an NSF CAREER Award, two Best Paper awards, and the UCSB Harold J. Plous Award (Junior Faculty of the Year).

I received a BS in Computer Science in 2001 from Cornell and a PhD in Computer Science from UNC Chapel Hill in 2006 under the supervision of Ming Lin. I was a Post-Doctoral Fellow at IBM TJ Watson Research Center in 2007, and a Post-Doctoral Associate at Cornell University from 2008-2009 under the supervision of Doug James.
From 2009-2011, I was an Assistant Professor in CS at the University of Saskatchewan. I was at UCSB from 2011-2017.
I interned for Rhythm and Hues Studios in 2000 and 2001, and my work was used for the Sorting Hat in Harry Potter and the Sorcerer's Stone.

The Laplacian Eigenfunction method for fluid simulation, which we refer to as Eigenfluids, introduced an elegant new way to capture intricate fluid flows with near-zero viscosity. However, the approach does not scale well, as the memory cost grows prohibitively with the number of eigenfunctions. The method also lacks generality, because the dynamics are constrained to a closed box with Dirichlet boundaries, while open, Neumann boundaries are also needed in most practical scenarios. To address these limitations, we present a set of analytic eigenfunctions that supports uniform Neumann and Dirichlet conditions along each domain boundary, and show that by carefully applying the discrete sine and cosine transforms, the storage costs of the eigenfunctions can be made completely negligible. The resulting algorithm is both faster and more memory-efficient than previous approaches, and able to achieve lower viscosities than similar pseudo-spectral methods. We are able to surpass the scalability of the original Laplacian Eigenfunction approach by over two orders of magnitude when simulating rectangular domains. Finally, we show that the formulation allows forward scattering to be directed in a way that is not possible with any other method.

Generating realistic fluid simulations remains computationally expensive, and animators can expend enormous effort trying to achieve a desired motion. To reduce such costs, several methods have been developed in which high-resolution turbulence is synthesized as a post process.
Since global motion can then be obtained using a fast, low-resolution simulation, less effort is needed to create a realistic animation with the desired behavior.
While much research has focused on accelerating the low-resolution simulation, the problem controlling the behavior of the turbulent, high-resolution motion has received little attention.
In this paper, we show that style transfer methods from image editing can be adapted to transfer the turbulent style of an existing fluid simulation onto a new one.
We do this by extending example-based image synthesis methods to handle velocity fields using a combination of patch-based and optimization-based texture synthesis. Importantly, this approach allows us to incorporate the incompressibility condition.
Using our method, a user can easily and intuitively create high-resolution fluid animations that have a desired turbulent motion.

Simulation artists frequently work with characters that self-intersect. When these characters are sent as inputs to a cloth simulator, the results can often contain terrible artifacts that must be addressed by tediously sculpting either the input characters or the output cloth. In this talk, we apply volume simulation to character meshes and remove self-intersections before they are sent to the cloth simulator. The technique has successfully dealt with very challenging animation scenarios in a production setting, and was applied to all the characters in the short film Bao.

Robustly simulating the dynamics of skin sliding over a character's body is an ongoing challenge. Skin can become non-physically "snagged" in curved or creased regions, such as armpits, and create unusable results. These problems usually arise when it becomes ambiguous which kinematic surface the skin should be sliding along. We have found that many of these problems can be addressed by performing 2D ray-tracing over the surface of the mesh. The approach is fast and robust, and has been used successfully in Incredibles 2.

Non-linear hyperelastic energies play a key role in capturing the fleshy appearance of virtual characters. Real-world, volume-preserving biological tissues have Poisson's ratios near 1/2, but numerical simulation within this regime is notoriously challenging. In order to robustly capture these visual characteristics, we present a novel version of Neo-Hookean elasticity. Our model maintains the fleshy appearance of the Neo-Hookean model, exhibits superior volume preservation, and is robust to extreme kinematic rotations and inversions. We obtain closed-form expressions for the eigenvalues and eigenvectors of all of the system's components, which allows us to directly project the Hessian to semi-positive-definiteness, and also leads to insights into the numerical behavior of the material. These findings also inform the design of more sophisticated hyperelastic models, which we explore by applying our analysis to Fung and Arruda-Boyce elasticity. We provide extensive comparisons against existing material models.

The intricate shapes and sounds that arise from vibrating Chladni plates are a well-known phenomenon. They are also quantitatively well understood, as the spatial patterns correspond to the eigenvectors of the underlying plate, and the audio frequencies arise from the plate's eigenvalues. We explore a generalization of the phenomenon by computing analogous quantities for a computational fluid dynamics simulation. Unlike the Chladni plate case, direct analytic expressions are not available, so we instead compute a set of "empirical" eigenvectors and eigenvalues. We find that these vectors form abstract, turbulent patterns in space. In another departure from the Chladni plate case, the eigenvalues no longer have a natural sonic mapping, so we construct a sonification that allows us to "listen" to the eigenvectors of the fluid. The united visual and sonic forms comprise a multimodal compositional palette that has great artistic potential.

We present a new method that achieves a two-way coupling between deformable solids and an incompressible fluid where the underlying geometric representation is entirely Eulerian. Using the recently developed Eulerian Solids approach [Levin et al. 2011], we are able to simulate multiple solids undergoing complex, frictional contact while simultaneously interacting with a fluid. The complexity of the scenarios we are able to simulate surpasses those that we have seen from any previous method. Eulerian Solids have previously been integrated using explicit schemes, but we develop an implicit scheme that allows large time steps to be taken. The incompressibility condition is satisfied in both the solid and the fluid, which has the added benefit of simplifying collision handling.

We propose a method to simulate the rich, scale-dependent dynamics of water waves. Our method preserves the dispersion properties of real waves, yet it supports interactions with obstacles and is computationally efficient. Fundamentally, it computes wave accelerations by way of applying a dispersion kernel as a spatially variant filter, which we are able to compute efficiently using two core technical contributions. First, we design novel, accurate, and compact pyramid kernels which compensate for low-frequency truncation errors. Second, we design a shadowed convolution operation that efficiently accounts for obstacle interactions by modulating the application of the dispersion kernel. We demonstrate a wide range of behaviors, which include capillary waves, gravity waves, and interactions with static and dynamic obstacles, all from within a single simulation.

Subspace fluid simulations, also known as reduced-order simulations, can be extremely fast, but also require basis matrices that consume an enormous amount of memory. Motivated by the extreme sparsity of Laplacian eigenfunctions in the frequency domain, we design a frequency-space codec that is capable of compressing basis matrices by up to an order of magnitude. However, if computed naively, decompression can be highly inefficient and dominate the running time, effectively negating the advantage of the subspace approach. We show how to significantly accelerate the decompressor by performing the key matrix-vector product in the sparse frequency domain. Subsequently, our codec only adds a factor of three or four to the overall runtime. The compression preserves the overall quality of the simulation, which we show in a variety of examples.

We present a method to increase the apparent resolution of particle-based liquid simulations. Our method first outputs a dense, temporally coherent, regularized point set from a coarse particle-based liquid simulation. We then apply a surface-only Lagrangian wave simulation to this high-resolution point set. We develop novel methods for seeding and simulating waves over surface points, and use them to generate high-resolution details. We avoid error-prone surface mesh processing, and robustly propagate waves without the need for explicit connectivity information. Our seeding strategy combines a robust curvature evaluation with multiple bands of seeding oscillators, injects waves with arbitrarily fine-scale structures, and properly handles obstacle boundaries. We generate detailed fluid surfaces from coarse simulations as an independent post-process that can be applied to most particle-based fluid solvers.

Subspace deformable body simulations can be very fast, but can behave unrealistically when behaviors outside the prescribed subspace, such as novel external collisions, are encountered. We address this limitation by presenting a fast, flexible new method that allows full space computation to be activated in the neighborhood of novel events while the rest of the body still computes in a subspace. We achieve this using a method we call subspace condensation, a variant on the classic static condensation precomputation. However, instead of a precomputation, we use the speed of subspace methods to perform the condensation at every frame. This approach allows the full space regions to be specified arbitrarily at runtime, and forms a natural two-way coupling with the subspace regions. While condensation is usually only applicable to linear materials, the speed of our technique enables its application to non-linear materials as well. We show the effectiveness of our approach by applying it to a variety of articulated character scenarios.

We present the first 3D algorithm capable of answering the question: what would a Mandelbrot-like set in the shape of a bunny look like? More concretely, can we find an iterated quaternion rational map whose potential field contains an isocontour with a desired shape? We show that it is possible to answer this question by casting it as a shape optimization that discovers novel, highly complex shapes. The problem can be written as an energy minimization, the optimization can be made practical by using an efficient method for gradient evaluation, and convergence can be accelerated by using a variety of multi-resolution strategies. The resulting shapes are not invariant under common operations such as translation, and instead undergo intricate, non-linear transformations.

We present an efficient new subspace method for simulating the self-contact of articulated deformable bodies, such as characters. Self-contact is highly structured in this setting, as the limited space of possible articulations produces a predictable set of coherent collisions. Subspace methods can leverage this coherence, and have been used in the past to accelerate the collision detection stage of contact simulation. We show that these methods can be used to accelerate the entire contact computation, and allow self-contact to be resolved without looking at all of the contact points. Our analysis of the problem yields a broader insight into the types of non-linearities that subspace methods can efficiently approximate, and leads us to design a pose-space cubature scheme. Our algorithm accelerates self-contact by up to an order of magnitude over other subspace simulations, and accelerates the overall simulation by two orders of magnitude over full-rank simulations. We demonstrate the simulation of high resolution (100K - 400K elements) meshes in self-contact at interactive rates (5.8 - 50 FPS).

We present a new subspace integration method that is capable of efficiently adding and subtracting dynamics from an existing high-resolution fluid simulation. We show how to analyze the results of an existing high-resolution simulation, discover an efficient reduced approximation, and use it to quickly "re-simulate" novel variations of the original dynamics. Prior subspace methods have had difficulty re-simulating the original input dynamics because they lack efficient means of handling semi-Lagrangian advection methods. We show that multi-dimensional cubature schemes can be applied to this and other advection methods, such as MacCormack advection. The remaining pressure and diffusion stages can be written as a single matrix-vector multiply, so as with previous subspace methods, no matrix inversion is needed at runtime. We additionally propose a novel importance sampling-based fitting algorithm that asymptotically accelerates the precomputation stage, and show that the Iterated Orthogonal Projection method can be used to elegantly incorporate moving internal boundaries into a subspace simulation. In addition to efficiently producing variations of the original input, our method can produce novel, abstract fluid motions that we have not seen from any other solver.

We propose a method of increasing the apparent spatial resolution of an existing liquid simulation. Previous approaches to this "up-resing" problem have focused on increasing the turbulence of the underlying velocity field. Motivated by measurements in the free surface turbulence literature, we observe that past certain frequencies, it is sufficient to perform a wave simulation directly on the liquid surface, and construct a reduced-dimensional surface-only simulation. We sidestep the considerable problem of generating a surface parameterization by employing an embedding technique known as the Closest Point Method (CPM) that operates directly on a 3D extension field. The CPM requires 3D operators, and we show that for surface operators with no natural 3D generalization, it is possible to construct a viable operator using the inverse Abel transform. We additionally propose a fast, frozen core closest point transform, and an advection method for the extension field that reduces smearing considerably. Finally, we propose two turbulence coupling methods that seed the high resolution wave simulation in visually expected regions.

Over the last decade, the special effects industry has embraced physics simulations as a highly useful tool for creating realistic scenes ranging from a small camp fire to the large scale destruction of whole cities. While fluid simulations are now widely used in the industry, it remains inherently difficult to control large scale simulations, and there is an constant struggle for increasing visual detail.

In this course, we will tackle these problems using turbulence methods. Turbulent detail is what makes typical fluid simulations look impressive, and the underlying physics motivate a powerful approach for control: they allow for an elegant split of large scale motion and small scale turbulent detail. This results in a two-stage work flow that is highly convenient for artists: first, a rough, and fast initial simulation is performed, which is then turned into a more detailed one by adding turbulent effects.

This course aims at giving an overview and providing practical guide to employing turbulence modeling techniques for fluid simulations in computer graphics. After reviewing the basics of fluid solvers, and the popular wavelet turbulence approach, we will present several powerful methods to capture advanced effects such as boundary layers, and turbulence with directional preferences. In addition, the difficulties of liquid simulations will be explained, and an approach for liquid turbulence that is based on wave dynamics will be presented.

We propose a domain-decomposition method to simulate articulated deformable characters entirely within a subspace framework. The method supports quasistatic and dynamic deformations, nonlinear kinematics and materials, and can achieve interactive time-stepping rates. To avoid artificial rigidity, or "locking," associated with coupling low-rank domain models together with hard constraints, we employ penalty-based coupling forces. The multi-domain subspace integrator can simulate deformations efficiently, and exploits efficient subspace-only evaluation of constraint forces between rotated domains using the so-called Fast Sandwich Transform (FST). Examples are presented for articulated characters with quasistatic and dynamic deformations, and interactive performance with hundreds of fully coupled modes. Using our method, we have observed speedups of between three and four orders of magnitude over full-rank, unreduced simulations.

Finite element simulations of nonlinear deformable models are computationally costly, routinely taking hours or days to compute the motion of detailed meshes. Dimensional model reduction can make simulations orders of magnitude faster, but is unsuitable for general deformable body simulations because it requires expensive precomputations, and it can suppress motion that lies outside the span of a pre-specified low-rank basis. We present an online model reduction method that does not have these limitations. In lieu of precomputation, we analyze the motion of the full model as the simulation progresses, incrementally building a reduced-order nonlinear model, and detecting when our reduced model is capable of performing the next timestep. For these subspace steps, full-model computation is "skipped" and replaced with a very fast (on the order of milliseconds) reduced order step. We present algorithms for both dynamic and quasistatic simulations, and a "throttle" parameter that allows a user to trade off between faster, approximate previews and slower, more conservative results. For detailed meshes undergoing low-rank motion, we have observed speedups of over an order of magnitude with our method.

We present a novel wavelet method for the simulation of fluids at high spatial resolution. The algorithm enables large- and small-scale detail to be edited separately, allowing high-resolution detail to be added as a post-processing step. Instead of solving the Navier-Stokes equations over a highly refined mesh, we use the wavelet decomposition of a low-resolution simulation to determine the location and energy characteristics of missing high-frequency components. We then synthesize these missing components using a novel incompressible turbulence function, and provide a method to maintain the temporal coherence of the resulting structures. There is no linear system to solve, so the method parallelizes trivially and requires only a few auxiliary arrays. The method guarantees that the new frequencies will not interfere with existing frequencies, allowing animators to set up a low resolution simulation quickly and later add details without changing the overall fluid motion.

We perform a detailed flop and bandwidth analysis of Jos Stam's Stable Fluids algorithm on the CPU, GPU, and Cell. In all three cases, we find that the algorithm is bandwidth bound, with the cores sitting idle up to 96% of the time. Knowing this, we propose two modifications to accelerate the algorithm. First, a Mehrstellen discretization for the pressure solver which reduces the running time of the solver by a third. Second, a static caching scheme that eliminates roughly 99% of the random lookups in the advection stage. We observe a 2x speedup in the advection stage using this scheme. Both modifications apply equally well to all three architectures.

Recent efforts to visually capture the phenomena of boiling have proposed monolithic approaches that extend the basic techniques underlying existing fluid solvers. In this work, we show that if we instead treat boiling as a separate computational module to be loosely coupled to an existing solver, a very easy to implement, highly efficient algorithm can be designed that produces excellent visual results, even on coarse (643) grids. The algorithm is also highly SIMD-amenable, allowing the boiling computation to be farmed out to a GPU or Playstation 3 Cell processor. Our algorithm takes less than 100 lines of commented, readable C++, and can be integrated into an existing particle level set fluid solver with virtually no modifications. A serial implementation consumes between 3-5% of the overall running time, and a preliminary SIMD implementation shows that a 643 simulation runs at 130 FPS, making the computational cost of the module totally negligible.

We present a technique for synthesizing spatially and temporally varying textures on continuous flows using image or video input, guided by the physical characteristics of the fluid stream itself. This approach enables the generation of realistic textures on the fluid that correspond to the local flow behavior, creating the appearance of complex surface effects, such as foam and small bubbles. Our technique requires only a simple specification of texture behavior, and automatically generates and tracks the features and texture over time in a temporally coherent manner. Based on this framework, we also introduce a technique to perform feature-guided video synthesis. We demonstrate our algorithm on several simulated and recorded natural phenomena, including river streams and lava flows. We also show how our methodology can be extended beyond realistic appearance synthesis to more general scenarios, such as temperature-guided synthesis of complex surface phenomena over a liquid during boiling.

Turing first theorized that many biological patterns arise through the processes of reaction and diffusion. Subsequently, reaction-diffusion systems have been studied in many fields, including computer graphics. We first show that for visual simulation purposes, reaction-diffusion equations can be made unconditionally stable using a variety of straightforward methods. Second, we propose an anisotropy embedding that significantly expands the space of possible patterns that can be generated. Third, we show that by adding an advection term, the simulation can be coupled to a fluid simulation to produce visually appealing flows. Fourth, we couple fast marching methods to our anisotropy embedding to create a painting interface to the simulation. Unconditional stability to maintained throughout, and our system runs at interactive rates. Finally, we show that on the Cell processor, it is possible to implement reaction-diffusion on top of an existing fluid solver with no significant performance impact.

We present a fast method for simulating, animating, and rendering lightning using adaptive grids. The "dielectric breakdown model" is an elegant algorithm for electrical pattern formation that we extend to enable animation of lightning. The simulation can be slow, particularly in 3D, because it involves solving a large Poisson problem. Losasso et al. recently proposed an octree data structure for simulating water and smoke, and we show that this discretization can be applied to the problem of lightning simulation as well. However, implementing the incomplete Cholesky conjugate gradient (ICCG) solver for this problem can be daunting, so we provide an extensive discussion of implementation issues. ICCG solvers can usually be accelerated using "Eisenstat's trick," but the trick cannot be directly applied to the adaptive case. Fortunately, we show that an "almost incomplete Cholesky" factorization can be computed so that Eisenstat's trick can still be used. We then present a fast rendering method based on convolution that is competitive with Monte Carlo ray tracing but orders of magnitude faster, and we also show how to further improve the visual results using jittering.

Laplacian instability is the physical mechanism that drives pattern formation in many disparate natural phenomena. However, current algorithms for simulating this instability are impractically slow and memory intensive. We present a new algorithm that is over three orders of magnitude faster than previous methods and decreases memory use by two orders of magnitude. Our algorithm is based on the dielectric breakdown model from physics, but is faster, more intuitive, easier to implement, and simpler to control. We demonstrate the ability of our algorithm to simulate various natural phenomena and compare its performance with previous techniques.

We present a novel technique for synthesizing textures over dynamically changing fluid surfaces. We use both image textures as well as bump maps as example inputs. Image textures can enhance rendering of the fluid by imparting novel realistic appearance to it, whereas bump maps enable the generation of complex micro-structures on the surface of the fluid that may be very difficult to synthesize using simulation. To generate temporally coherent textures over a fluid sequence, we transport texture information, i.e. color and local orientation, between fluid free surfaces from one time step to the next. This is accomplished by extending the texture information from the first fluid surface to the 3D fluid domain, advecting this information within the fluid domain along the fluid velocity field for one time step, and interpolating it back onto the second surface -- this operation, in part, uses a novel vector advection technique for transporting orientation vectors. We then refine the transported texture by performing texture synthesis over the second surface using our `surface texture optimization algorithm, which keeps the synthesized texture visually similar to the input texture and temporally coherent with the transported one. We demonstrate our novel algorithm for texture synthesis on dynamically evolving fluid surfaces in several challenging scenarios.

Large, 3D ice formations such as icicles exhibit a high degree of geometric and optical complexity. Modeling these features by hand can be a daunting task, so we present a novel physically-based algorithm for simulating this phenomenon. Solidification is usually posed as a so-called `Stefan problem', but the problem in its classic form is inappropriate for simulating the ice typically found in a winter scene. We instead use the `thin-film' variant of the Stefan problem to derive velocity equations for a level set simulation. However, due to the scales involved in the problem, even an adaptive grid level set solver is still insufficient to track the tip of an icicle. Therefore, we derive an analytical solution for the icicle tip and use it to correct the level set simulation. The results appear to be in agreement with experimental data. We also present a physically-based technique for modeling ripples along the ice surface that alleviates the need to explicitly track small-scale geometry. To our knowledge, our approach is the most complete model available, and produces complex visual phenomena that no previous method has been able to capture.

We present a physically-based method for animating and rendering lightning and other electric arcs. For the simulation, we present the dielectric breakdown model, an elegant formulation of electrical pattern formation. We then extend the model to animate a sustained, "dancing" electrical arc, by using a simplified Helmholtz equation for propagating electromagnetic waves. For rendering, we use a convolution kernel to produce results competitive with Monte Carlo ray tracing. Lastly, we present user parameters for manipulation of the simulation patterns.

We present a novel algorithm that simulates ice formation. Motivated by the physical process of ice growth, we develop a novel hybrid algorithm by synthesizing three techniques: diffusion limited aggregation, phase field methods, and stable fluid solvers. Each technique maps to one of the three stages of solidification. The visual realism of the resulting algorithm appears to surpass that of each technique alone, particularly in animations of freezing. In addition, we present a faster, simplified phase field method, as well as a unified parameterization that enables artistic manipulation of the simulation. We illustrate the results on arbitrary 3D surfaces.

The beautiful, branching structure of ice is one of the most striking visual phenomena of the winter landscape. Yet there is little study about modeling this effect in computer graphics. In this paper, we present a novel approach for visual simulation of ice growth. We use a numerical simulation technique from computational physics, the "phase field method," and modify it to allow aesthetic manipulation of ice crystal growth. We present acceleration techniques to achieve interactive simulation performance, as well as a novel geometric sharpening algorithm that removes some of the smoothing artifacts from the implicit representations. We have successfully applied this approach to generate ice crystal growth on 3D object surfaces in several scenes.

The Lévy Dragon is a well-known fractal introduced by Paul Lévy in 1938. It is a connected subset of the plane with interior (in fact it tiles the plane) but the interior is disconnected. Although the dragon has a fractal boundary of dimension 1.934007..., we show that each component of the interior has a polygonal boundary (with perhaps infinitely many edges) of finite length. There are infinitely many components, but we conjecture that they are all similar to one of sixteen different shapes. We show pictures of these shapes and some of the ways they interweave when two smaller dragons combine to make a larger dragon. We explain how we used the computer as a kind of "microscope" to reveal this structure. More pictures and programs are available on the web site http://www.mathlab.cornell.edu/~twk6/.

The geometric and optical complexity of ice has been a constant source of wonder and inspiration for scientists and artists. It is a defining seasonal characteristic, so modeling it convincingly is a crucial component of any synthetic winter scene. Like wind and fire, it is also considered elemental, so it has found considerable use as a dramatic tool in visual effects. However, its complex appearance makes it difficult for an artist to model by hand, so physically-based simulation methods are necessary. In this dissertation, I present several methods for visually simulating ice formation. A general description of ice formation has been known for over a hundred years and is referred to as the Stefan Problem. There is no known general solution to the Stefan Problem, but several numerical methods have successfully simulated many of its features. I will focus on three such methods in this dissertation: phase field methods, diffusion limited aggregation, and level set methods. Many different variants of the Stefan problem exist, and each presents unique challenges. Phase field methods excel at simulating the Stefan problem with surface tension anisotropy. Surface tension gives snowflakes their characteristic six arms, so phase field methods provide a way of simulating medium scale detail such as frost and snowflakes. However, phase field methods track the ice as an implicit surface, so it tends to smear away small-scale detail. In order to restore this detail, I present a hybrid method that combines phase fields with diffusion limited aggregation (DLA). DLA is a fractal growth algorithm that simulates the quasi-steady state, zero surface tension Stefan problem, and does not suffer from smearing problems. I demonstrate that combining these two algorithms can produce visual features that neither method could capture alone. Finally, I present a method of simulating icicle formation. Icicle formation corresponds to the thin-film, quasi-steady state Stefan problem, and neither phase fields nor DLA are directly applicable. I instead use level set methods, an alternate implicit front tracking strategy. I derive the necessary velocity equations for level set simulation, and also propose an efficient method of simulating ripple formation across the surface of the icicles.