  The goal of this paper is to create a fruitful bridge between the numerical
methods for approximating partial differential equations (PDEs) in fluid
dynamics and the (iterative) numerical methods for dealing with the resulting
large linear systems. Among the main objectives are the design of new efficient
iterative solvers and a rigorous analysis of their convergence speed. The link
we have in mind is either the structure or the hidden structure that the
involved coefficient matrices inherit, both from the continuous PDE and from
the approximation scheme: in turn, the resulting structure is used for deducing
spectral information, crucial for the conditioning and convergence analysis,
and for the design of more efficient solvers.
  As specific problem we consider the incompressible Navier-Stokes equations,
as numerical technique we consider a novel family of high order accurate
Discontinuous Galerkin methods on staggered meshes, and as tools we use the
theory of Toeplitz matrices generated by a function (in the most general block,
multi-level form) and the more recent theory of Generalized Locally Toeplitz
matrix-sequences. We arrive at a quite complete picture of the spectral
features of the underlying matrices and this information is employed for giving
a forecast of the convergence history of the conjugate gradient method,
together with a discussion on new more advanced techniques (involving
preconditioning, multigrid, multi-iterative solvers). Several numerical tests
are provided and critically illustrated in order to show the validity and the
potential of our analysis.
