Abstract: We present a highly efficient numerical solver for the Poisson equation on irregular voxelized domains supporting an arbitrary mix of Neumann and Dirichlet boundary conditions. Our approach employs a multigrid cycle as a preconditioner for the conjugate gradient method, which enables the use of a lightweight, purely geometric multigrid scheme while drastically improving convergence and robustness on irregular domains. Our method is designed for parallel execution on shared-memory platforms and poses modest requirements in terms of bandwidth and memory footprint. Our solver will accommodate as many as 7682 x 1152 voxels with a memory footprint less than 16 GB, while a full smoke simulation at this resolution fits in 32 GB of RAM. Our preconditioned conjugate gradient solver typically reduces the residual by one order of magnitude every 2 iterations, while each PCG iteration requires approximately 6.1 sec on a 16-core SMP at 7683 resolution. We demonstrate the efficacy of our method on animations of smoke flow past solid objects and free surface water animations using Poisson pressure projection at unprecedented resolutions.