Saturday, November 3, 2012

Analyzing Granular Flow Simulations in Mathematica

For the past two years, I have worked as a "data analyst" of sorts for the physics department at my university studying the results of numerical simulations of dense, gravity-driven, granular flows in a two-dimensional hopper, as in an hourglass. In this post, I talk a little bit about this research and share some of the Mathematica techniques that I have used to visualize numerical data.

The simulation that we use is event-driven and purely collisional. Event-driven means that the simulation "fast-forwards" in time to the next grain-grain or grain-boundary collision, calculates the new velocities (according to a partially inelastic interaction), and repeats this process continually. Purely collisional means that the grains do not form temporally-extended contacts. Because of the constant jostling about of grains, therefore, and because of the potential for inelastic collapsethe attempt of the simulation to resolve a cripplingly large number of collisions over the course of a short time intervalthe observation of "jamming" (i.e. sudden dynamical arrest) is extremely unlikely. Despite these odds, with careful sensitivity to several important parameters, we simulated an event strikingly close to a permanent jam:


You can watch the entire jamming transition in stunning HD here. I generated each frame from comma-separated-value (CSV) text file containing the following information for each of n=1000 particles:

particle 1: diameter | x position | y position | x velocity | y velocity 
particle n: diameter | x position | y position | x velocity | y velocity 

I reconstructed the spatio-temporal arrangement of the grains by inserting the x position, y position, and diameter of each grain (we use two different grain sizes: 1.0 and 1.2) into the following function: 

color [ diameter_ , xPos_ , yPos_ ] :=
   If [ diameter == 1.0 ,
      (* then *)
      Graphics [ { Gray , Disk [ { xPos , yPos } , diameter / 2 ] } ] ,
      (* else *)
      Graphics [ { Green , Disk [ xPos , yPos } , diameter / 2 ] } ] 

       ] ;

Looping this process for each of the thousands of files in the vicinity of the jam yields a collection of snapshots readily exportable to movie format with the simple command: 

Export [ " movie.avi " , listOfSnapshots , " FrameRate " → 60 ]

I incorporated the boundaries of the hopper into the visualization by calculating the relevant trigonometric relationships between geometrical parameters of the simulation such as top width, outlet width, taper angle, etc.

Things to Notice:

Grains near the center of the hopper travel faster in the direction of the driving force (i.e. gravity) than grains near the hopper walls. This is in marked contrast to, say, the homogeneous "plug flow" of water through an identical geometry. In granular flows, friction at the system boundaries produces inhomogeneity in the velocity field and induces the formation of "force chains" anchored at the hoppers walls that help support the weight of the granular column. This image from a Duke University experiment exemplifies this interesting phenomenon; the level of brightness (due to photoelasticity) indicates the magnitude of inter-grain stress:

More to come...

Sunday, October 7, 2012

Random Image Stitching

Panorama of my room created using the free open source, image stitching software Hugin. Photos taken with iPhone 4 and edited in GIMP.

Thursday, August 16, 2012

Buddhabrot: Deep Iteration

The Buddhabrot is an orbit density plot of points c for which the iteratively defined Mandelbrot sequence zn+1 = zn2 + c tends to infinity (z:= 0). It is a two-dimensional histogram that, roughly, represents the probability of a random orbit passing through a region(x,y) in the complex plane.

The Basic Technique

1. Select a random point c in the 4x4 square region centered at the origin of the complex plane.
2. Iterate c once according to the Mandelbrot equation. If the magnitude of z1 is less than two, compute z2. Repeat until either the absolute value of zn is no longer less than 2 or the iteration index is equal to a threshold iteration value maxDepth above which c is assumed (though not guaranteed) to remain bounded. The final iteration index corresponds to the orbit length of c
3. If the orbit length of c is less than maxDepth, recompute the orbit, this time mapping each point along the orbit path to an element of a res by res matrix, incrementing each "hit" matrix element (pixel) by +1
4. Repeat for many random points.
5. Convert this matrix to an image. I use the function ArrayPlot in Mathematica. 

This basic technique produces an image that, according to Wikipedia, resembles "classical depictions of Gautama Buddha."


1. Immediately identify and disregard initial points c that appear either in the cardioid (body) or in the main bulb (head) of the Mandelbrot set using the helpful equations on this page under the subheading "Optimizations." You can define more known M-set regions "by hand." (The irregular shape of the M-set, however, will limit further optimization using this technique to the realm of approximation.)
2. The Buddhabrot is symmetrical across the real axis. Mirror it for a speedup. I recommend avoiding this technique for deep iteration renderings, however, where asymmetrical orbits add to the aesthetic beauty of the image.
3. To avoid recomputing the orbits of good points, track and assign the orbit of every pointgood or badto a temporary res by res matrix. If the point turns out to be good, add the temporary matrix to the final image matrix. If not, throw the temporary matrix away. Warning: This method is a huge anti-optimization if one is seeking only high iteration points, since most temporary matrices will be throw-away matrices.

Deep Iteration

It is well known that filtering out the orbits of lower iteration points (very roughly ~10iterations) produces a Buddhabrot with beautiful and distinct orbits. The points that generate these sought-after orbits live very, very close to the border of the M-set. Finding these deep iteration points efficiently requires a clever sampling method.

Here is one possible method: 

1. Coarse-grain the 4x4 region centered at the origin of the complex plane into regularly spaced boxes; i.e. a grid.
2. For each point (x,y) on the grid, associate either +1 if orbitLength[(x,y)] is equal to maxBoxDepth or 0 if orbitLength[(x,y)] is less than maxBoxDepth. Remember: This algorithm only identifies general regions likely to yield high iteration points; there is no need to impose an ungodly high maximum iteration depth on the grid points themselves. 
3. Mark as a "good box" any set of four, rectangularly arranged grid points containing fewer than four but at least one +1. This amounts to identifying boxes with at least one corner, but not all corners, in the M-set. 

A finely-grained plane with a high maximum grid point iteration depth is ideal but computationally expensive. (Computation time increases with the square of the number of boxes per side.) The white boxes in the above plot come from a 100x100 grid with a maximum grid point iteration depth of 100. The boxes took only a few seconds to generate.

Using the above algorithmwith 3500 boxes per side and a maximum grid point iteration depth of two millionand over the course of about 24 hours, I found 297 points with iteration depths between 500k and one million. The left-hand figure is an array plot of the raw Buddhabrot data (the "hit" matrix) for a region in the belly of the Buddha; the right-hand figure is the same set of data third-root transformed.

The plot of the raw data is dark because, while a majority of pixels have similar, relatively low hit counts, a few particularly "attractive" pixels have very high hit counts. As a consequence, the coloring function stretches to accommodate these rare, high-hit-count pixels, leaving the low-hit-count, majority pixels black. The effect of the root transformation is to pull the hit counts closer together, but in a special way: The higher the hit count of a pixel, the further the transformation pulls it down. A quick plot of f(x)=x and f(x)=x1/3 shows that the distance between the two curves diverges with x. Here is the full (but scaled-down to 950x950) Buddhabrot:

The full-quality 5000x5000 Buddhabrot is available for download here. (No link yet.) Below is a zoom of the center of the Buddha. The brights have been colored orange and the mids have been colored magenta.