The thesis explores the interactions between cosmological gravitational waves (GWs) and large-scale magnetic fields. GWs are radiation produced by spacetime variations of the stress-energy tensor. Due to the weak coupling between gravity and the matter sector, GWs are a unique messenger from the early Universe before the Cosmic Microwave Background (CMB). Magnetic fields are observed across the Universe from the scales of planets and stars to galaxies and clusters, as well as the voids beyond the clusters. The present-day large-scale magnetic fields are believed to have evolved from primordial seed fields via magnetogenesis mechanisms active during the cosmic inflation and reheating epochs or the cosmological phase transitions occurring at the electroweak (EW) or quantum chromodynamic (QCD) scales in the early radiation-dominated (RD) era. The production of stochastic GW backgrounds (SGWBs) can be expected from the primordial electromagnetic (EM) fields or magnetohydrodynamic (MHD) turbulence around the time of EW and QCD phase transitions. The SGWBs then propagate through the pre-CMB Universe until the present day, carrying with them essential imprints of the corresponding sources as well as the underlying gravity theory at the early times. Since MHD turbulence is ubiquitously expected in the RD era, their induced SGWB spectrum is extensively studied.In one aspect of the thesis, we demonstrate that the MHD-GW system exhibits features of modified gravity (MG) in terms of the spectral slopes and amplitudes of the relic SGWB. We compute the spectra of GWs produced by MHD turbulence at the EW and QCD phase transitions, assuming massive gravity and scalar-tensor theories as two MG examples. Then we comparatively analyze these modified GW spectra with their counterparts in general relativity, and determine their qualitative and quantitative differences due to three effective MG parameters – graviton mass, GW friction, and GW speed. These spectral features are compared against the existing pulsar timing array (PTA) measurement in the nHz band as well as the expected GW detection sensitivities of upcoming instruments such as the Laser Interferometer Space Antenna (LISA) in the mHz band and Square Kilometer Array (SKA) as a PTA. The framework is general and can be applied to non-MHD sourced GW spectra. However, fixing MHD sources yields concrete constraints on the effective MG parameters.The other aspect of the thesis concerns the interaction between SGWBs from the early Universe and large-scale magnetic fields in the post-CMB Universe, which would convert a fraction of the gravitons in the SGWBs into photons of the same frequency via the inverse Gertsenshtein effect. The graviton-induced photons could then be a source of spectral distortions of the blackbody CMB.This, in turn, would allow us to constrain the GW amplitudes in the MHz-THz frequency regime. The high-frequency GWs (HFGWs) correspond to new physics in the early Universe such as beyond-Standard-Model mechanisms or sub-stellar mass primordial black holes. We place constraints on the HFGW energy densities by exploiting the existing direct observations of the radio sky, measurements of the 21-cm signal upper limits, the kinematic Sunyaev-Zeldovich observations, and assuming that graviton-induced photons saturate all of the reported radio excess over the CMB. We also forecast the potential of SKA and proposed future CMB surveys as novel HFGW detectors, and show that they will significantly tighten the current constraints and bring us a step closer to detecting HFGWs.