A domain decomposition spectral collocation method is developed for steady‐state, viscoelastic flow simulations through model porous media. The method is first tested on the one‐dimensional linear equation set resulting from the first‐order domain perturbation analysis of the flow of an upper convected Maxwell (UCM) fluid in an undulating channel. The splitting of the domain in the direction normal to the flow allows for an efficient local mesh refinement in the region of steep stress gradients and results in more accurate solutions than the fully pseudospectral implementation for the same degrees of freedom. The method was applied to the full two‐dimensional, nonlinear, viscoelastic equations with the same success. In addition to the undulating channel, solutions convergent with mesh refinement have been obtained for the flow of an UCM fluid through a square array of cylinders up toWe≊1.5. No increase in the flow resistance is found, for increasing elasticity. This observation is consistent with the findings in the undulating channel flow.