We are interested in parallel algorithms and software to enable accurate simulation of unsteady incompressible fluid flows in general three-dimensional geometries. The difficulty with this class of problems stems from a number of features of the governing Navier-Stokes equations. First, the highest order derivative is typically multiplied by a small number (the inverse of the Reynolds number). This singular perturbation results in thin boundary layers and gives rise to disparate length scales to be resolved by the numerical grid. Second, the equations are nonlinear, which leads to thin internal layers having unknown and possibly time varying positions. Finally, the incompressibility constraint must be satisfied at all times, implying global coupling of degrees-of-freedom at each time step.