This work has two parts. In the first, we treat a problem in evolving partial differential equations numerically. Typical methods for such problems are unstable if the time step is too big. The maximum allowed time step is limited to approximately the information propagation time between spatial grid points. Typical methods also use adaptive mesh refinement: for efficiency, the grid points are more finely spaced only in regions where the solution is rapidly varying. But then these regions have a smaller allowed time step. Using this small time step in regions where it is not needed is itself wasteful. While it would be better to be able to use the small time step only in regions where it is required, developing such local time-stepping methods can be difficult. We present a family of multistep integrators based on the Adams-Bashforth methods. These schemes can be constructed for arbitrary convergence order with arbitrary step size variation. The step size can differ between different subdomains of the system. It can also change with time within a given subdomain. The methods are linearly conservative, preserving a wide class of analytically constant quantities to numerical roundoff, even when numerical truncation error is significantly higher. These methods are intended for use in solving conservative PDEs in discontinuous Galerkin formulations, but are applicable to any system of ODEs. A numerical test demonstrates these properties and shows that significant speed improvements over the standard Adams-Bashforth schemes can be obtained. In the second part, we describe a new code, SpECTRE, for solving the GRMHD equations. This code uses the discontinuous Galerkin method and task-based parallelism to achieve scaling to exascale computing clusters. We have demonstrated that the code performs well on a variety of standard GRMHD test problems. We also show partial results from ongoing work evolving a relativistic disk surrounding a black hole.