In my last post, I wrote about a wave equation simulation that I spent a lot of time with as part of my ParallelAccelerator for Julia benchmarking effort. To illustrate the post, I made an animated GIF of the simulation running:
Over the course of the simulation, we can see circles emanating from the center of the plot, similar to how waves spread out in circles from the point where a drop of water hits a still pond. There are also waves coming from somewhere off to the left that interact with the waves coming from the center.
The simulation runs for 5000 time steps and then loops back to the beginning. There should be a frame of animation for each tenth time step, plus one initial frame — 501 frames in all. But something’s wrong — the animation looks choppy in places, as though some of the frames are missing. What gives?
I first noticed this frame-dropping behavior a long time ago, the first time I ran the code with visualization turned on. Because I was watching the visualization “live” as the code ran, though, I figured that my computer was just dropping frames because it was slow or busy with some other task. I thought the frames were there, just not being shown. It wasn’t until I made the above GIF and watched it that I suspected that something was wrong, because the GIF also looked like it was missing frames, even though it wasn’t running “live”. Hmmm.
I went back to look at the Julia code that did the plotting:
1 2 3 4
This code uses the Winston.jl 2D plotting package to plot an image for each tenth time step using the
imagesc function, then
displays that image on the screen. To make the above GIF, I’d also added an additional line of code that would save each image out to a file using Winston’s
savefig function, with sequentially numbered file names:
1 2 3 4 5
Since there’s a frame for every tenth time step, if we were on time step number 8, this code would save the image as
figure008.png, for instance.
Then I stitched together all the PNG files into an animation with a shell command using good ol’
The shell globbing conveniently takes care of putting the images in the correct order as frames in our GIF:
figure000.png, representing the initial state of the simulation, would be followed by
figure001.png and then
figure002.png, all the way up to
figure500.png. All this seemed to be working correctly; the frames weren’t being assembled in the wrong order. I still didn’t understand what was causing the jerkiness in the animation. But then I checked how many frames I actually had:
389? That was a lot fewer than the 501 I was supposed to have! I went back and looked at the plotting code again:
1 2 3 4 5
And then, finally, the problem jumped out at me. The line
if mod(t/dt, 10) == 0 is supposed to make sure that we only plot every tenth time step of the simulation.
t represents the amount of time that has passed since the beginning of the simulation, and
dt represents the amount of time that goes by with each time step. In my code, as in the original Octave code from which it was ported,
dt is set to 0.0001. Having
dt be a variable means that we can easily change the time resolution of the simulation by making
dt larger or smaller.
Why are we checking if
mod(t/dt, 10) is equal to zero? As the author of the original code explains: “
t/dt will return a natural number which represents the number of timesteps since the simulation has started. The modulo of this number and 10 will return the remainder of the division of these two numbers. Consequently comparing the modulo with 0 will only return true every tenth time the count variable
t was incremented […].”
The problem here is that, at least after porting to Julia from Octave,
t/dt can’t be counted on to return a natural number! Take, for instance, the third time step. At this time step,
t is 0.0003, and so
t/dt ought to be 0.0003 divided by 0.0001, which is 3. But because of floating-point inexactness,
t/dt comes out to something like 2.9999999999999996, and therefore
mod(t/dt, 10) also comes out to 2.9999999999999996.
The first dropped frame seems to occur on the 90th time step, when
t reaches 0.009. Then,
t/dt comes out to 89.99999999999999 instead of 90, and so
mod(t/dt, 10) comes out to 9.999999999999986 instead of 0. From that point on, frames are dropped here and there, whenever
t/dt works out to be something that isn’t an integer on a time step that happens to be a multiple of 10.
I tried logging the values of
mod(t/dt, 10) at each time step as the code ran. There were some stretches of time where
t/dt always came out integral:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
But during other stretches of time,
t/dt would be non-integral as often as not:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
Most interestingly of all, sometimes a pattern would emerge for a while — like here, for instance:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
And again, later on:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
During these stretches, we get an inexact answer for every fifth value of
t has a 2 or a 7 in its fourth decimal place, and not for any other value of
t. Other patterns show up elsewhere. Perhaps there’s something to be learned from studying these patterns. If you want to analyze the whole log of them, please be my guest!
Anyway, as is probably clear by now, we can fix the dropped-frames bug by rounding
t/dt up to the nearest integer:
1 2 3 4 5
With this fix, we have the correct number of frames:
convert command produces a lovely, smooth animation:
Ah! That’s better!
I updated my previous post to use the smooth animation instead of the choppy one. For comparison, here are the old and new versions, side by side. Since the smooth version has more than a hundred more frames than the choppy version, we can see that it also takes a bit longer to run:
Is there a lesson to be learned from all this? If the lesson of the previous post was that visualizing what your code is doing can be a helpful way to find bugs, perhaps the lesson here is to trust what the visualization shows you. In this case, our animation came out choppy for a reason. The other lesson here, I suppose, is to watch out for inexact fractions. In this case,
d/dt was already in the code we were porting from, but that’s no excuse!
Thanks to Cameron Finucane for a conversation that inspired this post.
Here’s something else that tripped me up at this stage: so,
-loopoption. I assumed that
-loop 1meant “loop” and
-loop 0meant “don’t loop”. But, no, in fact, the argument to
-loopis the number of times you want your GIF to loop – except in the special case of 0, which means “loop infinitely”. So
-loop 0had exactly the opposite of the behavior I expected! Oops. ↩