Rapidly oscillating Ap (roAp) stars exhibit an astrophysically interesting combination of strong, dipolar-like magnetic fields and high-overtone p-mode pulsations similar to the Sun. Recent time-resolved spectroscopy of these stars unravelled a complex picture of propagating magnetoacoustic pulsation waves, with amplitude and phase strongly changing as a function of atmospheric height. To interpret these observations and gain a new insight into the atmospheric dynamics of roAp stars we have carried out two-dimensional time-dependent, non-linear magnetohydrodynamical simulations of waves for a realistic atmospheric stratification of a cool Ap star. We explore a grid of simulations in a wide parameter space, treating oscillations of the velocity, magnetic field, and thermodynamic quantities in a self-consistent manner. Our simulations foster a new understanding of the influence of the atmosphere and the magnetic field on the propagation and reflection properties of magnetoacoustic waves, formation of node surfaces, and relative variation of different quantities. Our simulations reproduce all main features of the observed pulsational behavior of roAp stars. We show, for the first time, that the overall depth dependence of the pulsations in roAp atmospheres is strongly influenced by the density inversion at the photospheric base.