We present a novel method to simulate the propagation of seismic waves in realistic fluid-solid materials, coupled with dynamically evolving faults, in the self-gravitating prestressed Earth. A discontinuous Galerkin method is introduced, with a modified penalty numerical flux dealing with various boundary conditions, in particular with discontinuities. This numerical scheme allows general heterogeneity and anisotropy in the materials, by avoiding the diagonalization into polarized wave constituents such as in the approach based on solving elementwise Riemann problems, while maintains the numerical accuracy with mesh and polynomial refinements. We also include the interior slip boundary conditions for dynamic ruptures coupling with nonlinear friction laws, as an approach to simulate spontaneously cracking faults. We show the well-posedness for the system of particle motion coupled with gravitation field and its perturbation, by proving the coercivity of the bilinear operator, both in the continuous and discretized polynomial space, and therefore the convergence results. A multi-rate iterative scheme is proposed to address the challenging of solving the large implicit nonlinear system, and to allow different time steps for distinct physical processes in the overall coupling problem. We give rigorous proof for the well-posedness of mathematical model and moreover the stability of the numerical methods. Numerical experiments show the convergence as well as robustness in both well-established benchmark examples and realistic simulations.