We outline recent progress in understanding the accretion of plasma to rotating magnetized stars obtained from global axisymmetric (2D) and 3D magnetohydrodynamic (MHD) simulations in three main areas: (1.) Formation of jets from disk accretion onto rotating magnetized stars: From simulations where the viscosity and magnetic diffusivity within the disk are described by alpha models, we find long-lasting conical outflows/jets from the disk/magnetosphere boundary in both the case where the star is slowly rotating and where it is rapidly rotating (the “propeller regime”). Most of the mass flux in the outflows is in a hollow cone but inside this cone there is a low-density high-velocity magnetically dominated flow along the open polar field lines of the star. The outflows occur under conditions where the poloidal magnetic flux of the star is bunched up by the accretion disk near the disk/magnetosphere boundary. Recent simulations show that the conical outflows become well-collimated for axial distances of ≲ 20 times the inner disk radius. Exploratory 3D simulations show that conical winds are axisymmetric about the rotational axis (of the star and the disk), even when the dipole field of the star is significantly misaligned. (2.) Formation of intrinsically one-sided jets from disk accretion to rotating magnetized stars: There is strong observational evidence for an asymmetry between the approaching and receding jets from a number of young stars. We discuss the first MHD simulations of the formation asymmetric or one-sided jets arising from disk accretion to a rotating star with an asymmetric (dipole plus quadrupole) magnetic field. (3.) Global axisymmetric and 3D simulations of the magnetorotational instability (MRI) in disk accretion onto magnetized stars: In the axisymmetric simulations we observe cases where there is episodic or quasi-periodic burst of accretion similar to that observed in one X ray source. In 3D MHD simulations of accretion onto stars with tilted dipole fields using our Godunov-type code based on the “cubed sphere” grid we find that the density distribution is much less smooth than in the case of the laminar accretion flow described by α–viscosity. Instead, large turbulent cells dominate the flows and are strongly elongated in the azimuthal direction.