A high-order discretization consisting of a tensor product of the Fourier collocation and discontinuous Galerkin methods is presented for numerical modeling of magma dynamics. The physical model is an advection-reaction type system consisting of two hyperbolic equations and one elliptic equation. The high-order solution basis allows for accurate and efficient representation of compaction-dissolution waves that are predicted from linear theory. The discontinuous Galerkin method provides a robust and efficient solution to the eigenvalue problem formed by linear stability analysis of the physical system. New insights into the processes of melt generation and segregation, such as melt channel bifurcation, are revealed from two-dimensional time-dependent simulations.