Abstract: Most of the baryonic matter in the Universe exists in a state of turbulence. On Earth, turbulence causes atmospheric pressure fluctuations that make stars twinkle, and influences the chaotic movements of oceans and weather patterns. Turbulence is so prevalent on Earth that it is not surprising that it plays a fundamental role on almost all scales of the Universe. At the scale of hundreds of light years, turbulence prevails in the dusty, multiphase plasma environment found in the discs of galaxies, known as the interstellar medium (ISM). The ISM turbulence, filled with supersonic shocks and gas density fluctuations, and shaped by magnetic fields, differs significantly from terrestrial turbulence and lacks a robust theory. Throughout my PhD I embarked on an extensive exploration of ISM turbulence, focusing on three aspects: the turbulent dynamo, magnetised turbulent fluctuations, and the stationary turbulent energy cascade in supersonic, magnetised turbulence, relevant to the ISM. Using direct numerical simulations, I unveil universal long-term dynamics of the turbulent dynamo across different initial magnetic field conditions. Additionally, I craft detailed statistical models for various fluctuations in magnetohydrodynamic (MHD) turbulence, taking into account compressibility and anisotropy effects. These models, covering gas density, ion Alfvén velocity, and magnetic field fluctuations, have been extensively validated via numerical experiments. They have implications for star formation theories, cosmic ray transport, and contribute directly to current efforts to measure the interstellar magnetic field amplitude using dust polarisation measurements. Finally, I present the world's largest compressible MHD turbulence simulation, a simulation run for over 80 million core hours on approximately 140,000 compute cores. On larger scales, we encounter Burgers-like turbulence -- the turbulence of saw-tooth shocks, while on smaller scales, we find Boldyrev’s dynamic aligned turbulence. We also discern a power-law structure in the magnetic energy spectrum only at the highest magnetic Reynolds numbers, previously unexplained by any (even asymptotic) magnetic spectrum theories. This calculation, demonstrating such a dichotomy in the magnetised regime, bears significant implications for ISM turbulence and fundamental plasma and fluid sciences.