Obtaining accuracy along domain boundaries is often the primary goal of a numerical simulation. In this work, we show how the test space of discontinuous Galerkin (DG) methods can be optimized to achieve enhanced boundary accuracy for linear problems. Given some norm in which accuracy is desired, the optimal test functions render the numerical solution the best approximation to the true solution in that norm. For most norms, these test functions would be global in nature. However, for norms emphasizing boundary accuracy, they can be computed in an element-local manner and represent adjoint solutions for the interface fluxes. The resulting accuracy in these interface fluxes propagates globally, leading to convergence rates on the domain boundaries greater than 2p+1. Indeed, if the test functions and interface fluxes are well-resolved, exact boundary fluxes are obtained. Here, we demonstrate these ideas for several Computational Fluid Dynamics (CFD) simulations, including advection–diffusion and linearized Euler problems in both one and two dimensions. An extension of the theory to nonlinear problems will be presented in a future work.