We address in this thesis the current need to design new parallel algorithms and tools that ease the development of geodynamic modelling applications that are suited for today's and tomorrow's hardware. We present (1) the MATLAB HPC compiler HPC.m, which greatly simplifies the building of parallel high performance applications and (2) parallel algorithms for the 3D simulation of strongly nonlinear processes as mechanical and reactive porosity waves. To simulate mechanical porosity waves we employ a massively parallel algorithm that permits to resolve the deformation of fluid-filled viscoelastic porous media in 3D. The utilized mathematical model is based on Biot's poroelastic theory, extended to account for viscous deformation and plastic yielding. The modelling results exhibit the impact of decompaction weakening on the formation of three-dimensional solitary-wave-like moving porosity channels. To simulate reactive porosity waves we use a solver for 3D deformation of fluid-filled reactive viscous porous media. The Damköhler number (Da) of the simulations is varied in order to estimate the respective roles of viscous deformation (low Da) and reaction (high Da) on wave propagation. 3D waves are found to propagate independently of their source at constant speed by going through each other for all the investigated Da. Soliton-like wave propagation as a result of metamorphic reaction provides an efficient mechanism for fluid flow in the Earth's crust. We illustrate the great performance and versatility of HPC.m by deploying it to generate solvers for a variety of physics across multiple Earth Science disciplines. All solvers run close to hardware's peak performance and were shown to scale linearly on a institute cluster with 80 GPUs. Moreover, our nonlinear poroviscoelastic two-phase flow solver scales close to ideally on Piz Daint's 5000 GPUs at the Swiss National Supercomputing Centre.