Systems of hyperbolic conservation laws (HCLs) commonly arise as mathematical descriptors of the natural world, and are particularly ubiquitous in fluid dynamics. These laws appear as complicated and highly nonlinear partial differential equations describing the evolution of fundamental conserved quantities such as mass, momentum, and energy. Solving these equations analytically is entirely intractable for all but the simplest cases, and investigating problems with real world importance falls to numerical approaches more and more frequently. Most HCLs exhibit rich dynamics with complicated smooth flows and discontinuities coexisting, often with shocks arising frominitially smooth data. Designing numerical schemes that can efficiently and accurately represent smooth phenomena, while also remaining robust and reliable in the vicinity of shocks, is very challenging. Finite volume methods are one particularly useful approach to designing such methods as conservation is enforced discretely, and discontinuities can be represented quite naturally. An unfortunate drawback of these methods is that achieving high-order accuracy in multiple space dimensions is difficult. This dissertation overcomes these challenges by developing a kernel-based non-polynomial reconstruction scheme that is manifestly multidimensional. This scheme is first posed as a linear recovery problem in a reproducing kernel Hilbert space. This linear reconstruction method is then cast into a weighted essentially non-oscillatory (WENO) framework so that it may represent both smooth and discontinuous data. This scheme is then incorporated into solvers for the compressible Euler equations, compressible Navier-Stokes equations, and ideal magnetohydrodynamics (MHD) equations. In doing so, a novel set of variables that are more suited to multidimensional reconstruction, dubbed the linearized primitive variables, are introduced. Troubled cell indicators are developed that allow for a more accurate and efficient treatment of smooth solutions in an entirely automatic fashion. Positivity preserving limiters are also incorporated, and allow for the evolution of flows with extremely strong shocks. A highly parallel multi-GPU implementation is provided, and the proposed method is tested against a variety of stringent benchmark problems.