Optimization problems subject to PDE constraints form a mathematical tool that can be applied to a wide range of scientific processes, including fluid flow control, medical imaging, option pricing, biological and chemical processes, and electromagnetic inverse problems, to name a few. It is necessary to obtain accurate solutions to such problems within a reasonable computation time, in particular for time-dependent structures, for which the "all-at-once" solution can lead to extremely large linear systems. In this talk we consider iterative methods, in particular Krylov subspace methods, to solve systems of this form arising from fluid flow control problems, and we accelerate these methods using fast and robust preconditioning strategies. In particular, we are interested in devising periodic-in-time approximations of the huge-scale linear systems, to enable the solution of the problems via fast Fourier transforms and block diagonal solves. This leads to parallel-in-time preconditioned iterative methods. The system as a whole is solved using the flexible GMRES algorithm, with inner solvers involving the application of GMRES or MINRES. The diagonal solves are implemented using a variety of strategies, including optimal block diagonal preconditioners (for Stokes problems), Uzawa iteration, Chebyshev semi-iteration and multigrid methods for individual matrices, as well as recently-devised block commutator arguments for matrices discretizing "divergence x matrix-inverse x gradient" terms. We illustrate the benefit of our new approach through a range of numerical experiments for Stokes and Oseen problems, and in particular we observe convincing strong and weak scaling results.