In this paper, we study a two level method based on Newton's iteration for the nonlinear system arising from the Galerkin finite element approximation to the equations of motion described by the Kelvin-Voigt viscoelastic fluid flow model. The two-grid algorithm is based on three steps. In the first step, the nonlinear system is solved on a coarse mesh 𝒯 H $\mathcal {T}_{H}$ to obtain an approximate solution u H . In the second step, the nonlinear system is linearized around u H based on Newton's iteration and the linear system is solved on a finer mesh 𝒯 h $\mathcal {T}_{h}$ . Finally, in the third step, a correction to the results obtained in the second step is achieved by solving a linear problem with a different right hand side on 𝒯 h $\mathcal {T}_{h}$ . Optimal error estimates in L ? (L 2)-norm, when h = 𝒪 ( H 2 ? ? ) $h=\mathcal {O} (H^{2-\delta })$ and in L ? (1)-norm, when h = 𝒪 ( H 5 ? 2 ? ) $h=\mathcal {O}(H^{5-2\delta })$ for velocity and in L ? (L 2)-norm, when h = 𝒪 ( H 5 ? 2 ? ) $h=\mathcal {O}(H^{5-2\delta })$ for pressure are established, where ? > 0 is arbitrarily small for two dimensions and ? = 1 2 $\delta =\frac {1}{2}$ for three dimensions.