The hybridized discontinuous Galerkin methods (HDG) introduced a decade ago is a promising candidate for high-order spatial discretization combined with implicit/implicit-explicit time stepping. Roughly speaking, HDG methods combines the advantages of both discontinuous Galerkin (DG) methods and hybridized methods. In particular, it enjoys the benefits of equal order spaces, upwinding and ability to handle large gradients of DG methods as well as the smaller globally coupled linear system, adaptivity, and multinumeric capabilities of hybridized methods. However, the main bottleneck in HDG methods, limiting its use to small to moderate sized problems, is the lack of scalable linear solvers. In this thesis we develop fast and scalable solvers for HDG methods consisting of domain decomposition, multigrid and multilevel solvers/preconditioners with an ultimate focus on simulating large scale problems in fluid dynamics and magnetohydrodynamics (MHD). First, we propose a domain decomposition based solver namely iterative HDG for partial differential equations (PDEs). It is a fixed point iterative scheme, with each iteration consisting only of element-by-element and face-by-face embarrassingly parallel solves. Using energy analysis we prove the convergence of the schemes for scalar and system of hyperbolic PDEs and verify the results numerically. We then propose a novel geometric multigrid approach for HDG methods based on fine scale Dirichlet-to-Neumann maps. The algorithm combines the robustness of algebraic multigrid methods due to operator dependent intergrid transfer operators and at the same time has fixed coarse grid construction costs due to its geometric nature. For diffusion dominated PDEs such as the Poisson and the Stokes equations the algorithm gives almost perfect hp--scalability. Next, we propose a multilevel algorithm by combining the concepts of nested dissection, a fill-in reducing ordering strategy, variational structure and high-order properties of HDG, and domain decomposition. Thanks to its root in direct solver strategy the performance of the solver is almost independent of the nature of the PDEs and mostly depends on the smoothness of the solution. We demonstrate this numerically with several prototypical PDEs. Finally, we propose a block preconditioning strategy for HDG applied to incompressible visco-resistive MHD. We use a least squares commutator approximation for the inverse of the Schur complement and algebraic multigrid or the multilevel preconditioner for the approximate inverse of the nodal block. With several 2D and 3D transient examples we demonstrate the robustness and parallel scalability of the block preconditioner