Variability is a generic feature of accretion onto black holes. In both X-ray binaries and active galactic nuclei, variability is observed on nearly all accessible timescales and across the entire electromagnetic spectrum. On different timescales and at different wavelengths it has unique signatures that can be used to characterize the accretion processes generating the emission and probe the accretion disks, which would otherwise be impossible. Despite having been observed for over fifty years, interpreting this variability is difficult. Simple phenomenological models have been used to explain the behaviors and geometries of the observed accretion disk, but they have yet to be rigorously tested in a full magnetohydrodynamic framework. In this dissertation we use high-resolution numerical models to investigate: (1) ``propagating fluctuations" in mass accretion rate that give rise to the nonlinear signatures of accretion on viscous timescales, (2) the dynamics of truncated accretion disks which are invoked to explain the spectral variation of outbursting X-ray binaries and the bifurcation of AGN accretion states, and (3) the large-scale magnetic dynamo behavior in thick and thin accretion disks.

We find that the structured variability readily seen in the light curves from accreting black holes (i.e. log-normal flux distributions, linear relations between the RMS and the flux, and radial coherence) quickly and naturally grows from the MRI-driven turbulence and that these properties translate into photometric variability. For the first time, we identify the large-scale magnetic dynamo as the source of the low-frequency modulations of the disk stress that cause this structure. We introduce a bistable cooling law into hydrodynamic and magnetohydrodynamic simulations to study the manifestation of a truncated accretion disk in each regime. We find that rather than a truncation edge, the transition is better described by a truncation zone" when the angular momentum transport and heating is governed by MRI-driven turbulence instead of a true viscosity. Additionally, we find that the hot gas in the simulation buoyantly rises in a gentle outflow and eventually fills the entire volume, instead of simply being confined to the innermost region. The outflow interacts with the disk body and enhances the magnetic stresses, which could produce stronger quasiperiodic variability. Finally, we conduct an investigation of the large-scale magnetic dynamo using a suite of four global magnetohydrodynamic disk simulations with scaleheight ratios of $h/r=\{0.05, 0.1, 0.2, 0.4\}$. Most notably, the organization that is prevalent in accretion disk simulations and described as a butterfly pattern" does not occur when $h/r \ge 0.2$, despite the dynamo action still operating efficiently.