Plasmas are composed of a large number of charged and neutral particles that interact through electromagnetic forces and collisions over a wide range of temporal and spatial scales to produce collective dynamics. The discrete nature of the plasma particles, their large number, and the wide range of interaction scales make a many-body treatment computationally prohibitive. A statistical treatment evolves the probability density function for each particle species in phase space (x,v) and produces the more manageable but six-dimensional kinetic model, which is typically considered to provide the highest physical fidelity. Solutions to the kinetic model are possible by phase space discretizations using high-order finite element methods that exploit GPU-based HPC systems. Solutions using the discontinuous Galerkin method will be described, where the governing equations are expressed in a balance law form, an approximate Riemann solver is developed for evaluating inter-element numerical fluxes, and Runge-Kutta methods perform the temporal advance. Tailored initializations using kinetic eigenfunctions are shown to access saturated physics regimes with computational efficiency. Assumptions related to thermodynamic relaxation reduce the kinetic model through velocity moments to retain more limited information about the distribution function at each position (x). The moment models, e.g. 5N, 10N, and 13N, are three dimensional and thereby offer substantial computational acceleration compared to the kinetic model. The validity of the assumptions and the resulting moment models depends on local plasma parameters such as collisionality, charge neutrality, and magnetization. Such considerations are necessary when selecting the appropriate model for simulations. Expressing the hierarchy of continuum plasma models in a consistent formulation facilitates hybridization. The computational methods are applied to the GEM collisionless magnetic reconnection problem, 2D electrostatic and electromagnetic kinetic instabilities, and the magnetized Kelvin-Helmholtz instability. Spatially localized deviations from thermodynamic equilibrium are investigated, as well as the resulting agreement and disagreement of global solution metrics between plasma models of differing fidelity. The results suggest opportunities to hybridize by applying the simplest, locally valid plasma model and developing strategies to couple the models across subdomain boundaries.