A time-marching pressure-correction/Taylor-Galerkin finite element algorithm is presented to accommodate low Mach number compressible and incompressible viscoelastic liquid flows. The algorithm is based on an operator splitting constructive process that discloses three fractional stages. For the compressible regime, a piecewise-constant density interpolation with gradient recovery is employed, for which the background theory and consistency of approach are discussed. The scheme is applied to contraction flows for Oldroyd model fluids, covering entry-exit flows and high pressure-drop situations. Stability and performance characteristics of the new algorithmic implementation are highlighted. Solutions are provided for a range of compressible settings, tending to the incompressible limit at vanishing Mach number. (C) 2004 Elsevier B. V. All rights reserved.