We consider techniques that enable large-scale gradient-based shape optimization of wave-guiding devices in the context of three-dimensional time-domain simulations. The approach relies on a memory efficient boundary representation of the shape gradient together with primal and adjoint solvers semiautomatically generated by the FEniCS framework. The hyperbolic character of the governing linear wave equation, written as a first-order system, is exploited through systematic use of the characteristic decomposition both to define the objective function and to obtain stable numerical fluxes in the discontinuous Galerkin spatial discretization. The methodology is successfully used to optimize the shape of a midrange acoustic horn, described by 1,762 design variables, for maximum transmission efficiency, where the parallel computations involve a total of $3.5\times10^9$ unknowns.