We implement a high-order finite-element application, which performs the numerical simulation of seismic wave propagation resulting for instance from earthquakes at the scale of a continent, from active seismic acquisition experiments in the oil industry, or from ocean acoustics problems, on a large GPU cluster. Mesh coloring enables an efficient accumulation of degrees of freedom in the assembly process over an unstructured mesh. We use non-blocking MPI and show that computations and communications over the network and between the CPUs and the GPUs are almost fully overlapped. The MPI + GPU solver exhibits good scalability and achieves significant speedup over a carefully tuned equivalent CPU code. In terms of applications, the challenge is to solve large forward and adjoint/inverse elastic wave propagation problems at very high resolution, i.e., at very high frequency, keeping in mind that solving such inverse tomography problems iteratively implies running thousands of simulations for a given problem.