Nice writeup. One thing worth adding to the limitations: without vorticity confinement, the Gauss-Seidel projection step quietly dissipates the small-scale curl that makes smoke look like smoke.

The 2001 Fedkiw/Stam/Jensen "Visual Simulation of Smoke" paper added it back as a correction force for exactly this reason. At N=16 it doesn't matter much because the grid itself can't represent fine vortices, but the moment you crank N up the missing confinement becomes visible.