"If a picture is worth a thousand words, then a movie is worth a thousand pictures, a simulation is worth a thousand movies, and a framework is worth a thousand simulations."
Introduction
As part of my final project for a computational astrophysics course at the University of Arizona, I developed a 2D model of impact cratering. Typically, simulations of impact crater processes are conducted using fluid dynamics equations that capture shock waves that travel through the medium or by using advanced particle codes that . I wanted to see if a very simple 2D model could reproduce some of the results of these more sophisticated approaches. (I also wanted to try to implement several optimization techniques I thought about. More on that later.)
Animations
I experimented with a number of parameters including the ratio of the densities of the impactor and the ejecta, packing density (balls per unit area), impactor speed, ball size, impactor size, and more. Because I posted videos of these animations online, I have been asked a number of questions about them from people around the world. It is on my to-do list to make the code publicly available.
To make the animations, I tell the simulation how often to take a picture. The simulation then makes sure to stop to take a picture (saving it to a file) before continuing to iterate. With careful control of the timestepping — the simulation now has to think about stopping to handle a collision and stopping to take a picture — this ensures a smooth animation for the viewer. At the end, I combine the hundreds or thousands of images together in another program to make the video.
My Model
In my model, I treat the ground as a "ball pit" containing many small particles that are initially at rest. I treat the impactor as either one large ball or as a clump of many smaller balls, typically with a higher mass density than the ground. The impactor is given an initial velocity such that it will hit the ground. I treat collisions between particles as "hard shell", elastic collisions, where "hard shell" means that they do not crush or crumble during an impact, and "elastic" means that energy is perfectly conserved in every collision. Because these collisions happen on very short timescales, I ignore gravity.
I color the balls differently to differentiate the impactor material and the ground material. I added a "tracing" feature to some fraction of the ground balls to keep track of where they move on the screen. Unlike most of my Java simulations, this does not run fast enough to watch in real time unless the ball count is low (say, less than 5000). Instead, I save pictures and make a movie after the simulation runs.
In the videos that follow, you can see that the particles move downward (left on the screen) and outward (up and down on the screen) from the impact site, creating an underground shockwave that propagates at the sound speed. Particles near the surface generally move outward at a roughly 45–55° angle away from the surface. This is in excellent agreement with the predictions of the Maxwell Z-model (Maxwell, 1977), the most influential theoretical model of impact crater streamlines.
Algorithms
Because particle interactions are treated as hard-shell, elastic collisions, I can freeze time and calculate which particles will collide with which other particles and how much time will elapse until they collide. To do this, I ask each ball to calculate whether it will collide with each other ball, based on their positions and velocities, and assuming that they're the only two balls in existence. From there, I look at all possible collisions to figure out when the next collision will happen and which two balls are involved in that collision. I then fast-forward time to the next collision where the two balls are touching, handle the collision by updating their velocities (they bounce off each other), and then repeat the process for the next collision.
Asking each particle to calculate if and when it will collide with each other particle takes a very long time. If I have $N$ balls, I will need to compute $N \times (N-1)$ collisions, which is an $\mathcal{O}(N^2)$ operation. Because of this, I use a number of optimizations to speed up the overall calculations. Generally, this means being clever about not calculating things multiple times. I detail these tricks below.
Optimization #0: Avoid Redundancy
The easiest way to avoid calculating unnecessary collisions is to only calculate them once. If I ask Ball A whether it collides with Ball B, I do not need to later ask Ball B whether it collides with Ball A — I've already done that calculation. This reduces the number of collisions I need to calculate by half. (Unfortunately, looking back at my code, this step is not done!)
Optimization #1: Collision Queue
A common method for avoiding repeating calculations is to store the results of a calculation for later. To this end, I first directed each particle to come up with a list of all possible collisions, and then combined them into a master list (or rather, a queue) of all collisions that could happen, sorted by how soon they would happen. I call this the Collision Queue.
When two balls collide, they change their trajectories. This means I need to re-calculate their collision list because all of their potential collisions will no longer occur. So, I remove all possible collisions in the Collision Queue, then ask each of the two balls for an updated list of new possible collisions, and add these to the master Collision Queue. All the potential collisions between other balls are not affected and may still happen in the future. I therefore leave their collisions in the Collision Queue until they change course.
Optimization #2: Only Look for Next Collision
Because a ball changes its trajectory when it collides, I only need to ask each ball for its next collision — all subsequent collisions will be recalculated anyway. The main purpose of this optimization is to reduce the number of potential collisions in the master Collision Queue. A smaller Collision Queue makes it faster to insert new potential collisions or remove collisions that are no longer applicable (e.g. because one of the particles involved in that collision hit something else first).
Optimization #3: Grid Cells
It is very time-consuming for a ball to recalculate its list of collisions with all other balls. Even if only two balls have to do this with all other balls, this is still a lot of calculations for simulations with tens of thousands of balls. With so many balls, one might ask whether it's really necessary to for a ball on one side of the screen to check for a collision with a ball on the other side of the screen. After all, with so many balls in the way, shouldn't we expect them to collide with other balls first?
To reduce the number of calculations, I divide the entire simulation area into a 2D grid of rectangular cells of a carefully-chosen size. Particles in each cell only look for collisions with particles in their own cell and in neighboring cells. From the standpoint of object-oriented design, each grid cell keeps a list of which balls are inside of it. When a ball leaves one cell in the course of the simulation, ownership of that ball transfers to the new cell. At the cost of some overhead, this significantly cuts down the number of potential collisions each particle has to check for. The smaller the cells are, the more potential collisions can be ignored. Typically, I set the minimum cell size to be twice the meteor radius. Ignoring empty "above-ground" cells, most cells have only a couple particles. This was perhaps the single most effective optimization.
Future Work
I would like to try this again using MPI on a small cluster or using a GPU for acceleration. I would also really love to add some setup options for this model so others could play around with it and make their own movies with it. I would also like to make it a standalone app so teachers could use it in the classroom.