The Q^k spectral element method has been a popular high order method for solving second order PDEs for nearly four decades, obtained by continuous finite element method with tenor product polynomial of degree k basis and with (k+1)-point Gauss-Lobatto quadrature. In this talk, I will present some recent results of this classical scheme, including its accuracy, monotonicity (stability), and examples of using monotonicity to construct high order accurate bound (or positivity) preserving schemes in various applications including the Allen-Cahn equation, Keller-Segel equation for chemotaxis, and especially compressible Navier-Stokes equations for high speed flows. A simple and efficient implementation in MABLAB 2023 for GPU acceleration of Poisson solvers will also be presented.