I have fixed a number of bugs and now have a solver based on BFECC advection using min-max limited cubic interpolation for non-uniform and often highly anisotropic meshes for the velocity and pressure solve, with high-resolution uniform density fields for the density field. The results are fairly impressive:
It's still not perfect, but given the cell size and anisotropy, I think it does extremely well. Although there are definite artifacts, the payoff is in the memory usage and runtime. The whole simulation is 5 seconds in real time and takes approximately 22 seconds per output frame, meaning my 5 second simulation completes in under an hour. These simulations used to take on the order of 4-5 hours.
The results look pretty good. I think the grading is too steep to get really nice results, but it's an excellent proof-of-concept.