DescriptionMany scientific domains require the solution of high dimensional PDEs. Traditional grid- or mesh-based methods for solving such systems in a noise-free manner quickly become intractable due to the scaling of the degrees of freedom going as O(N^d) sometimes called "the curse of dimensionality." We are developing an arbitrarily high-order discontinuous-Galerkin finite-element solver that leverages an adaptive sparse-grid discretization whose degrees of freedom scale as O(N*log2 N^D-1). This method and its subsequent reduction in the required resources is being applied to several PDEs including time-domain Maxwell's equations (3D), the Vlasov equation (in up to 6D) and a Fokker-Planck-like problem in ongoing related efforts. Here we present our implementation which is designed to run on multiple accelerated architectures, including distributed systems. Our implementation takes advantage of a system matrix decomposed as the Kronecker product of many smaller matrices which is implemented as batched operations.