Fractional calculus and fractional-order modeling provide effective tools for modeling and simulation of anomalous diffusion with power-law scalings. In complex multi-fractal anomalous transport phenomena, distributed-order partial differential equations appear as tractable mathematical models, where the underlying derivative orders are distributed over a range of values, hence taking into account a wide range of multi-physics from ultraslow-to-standard-to-superdiffusion/wave dynamics. We develop a unified, fast, and stable Petrov-Galerkin spectral method for such models by employing Jacobi poly-fractonomials and Legendre polynomials as temporal and spatial basis/test functions, respectively. By defining the proper underlying distributed Sobolev spaces and their equivalent norms, we rigorously prove the well-posedness of the weak formulation, and thereby, we carry out the corresponding stability and error analysis. We finally provide several numerical simulations to study the performance and convergence of proposed scheme.