Simulation and analysis of stiff PDE requires robust and algorithmically scalable solvers. Multilevel methods are required to provide this algorithmic scalability, but appear in various forms due to additional problem structure, physical regime, coupling strengths, discretizations, problem size, and hardware. I will discuss software components developed in PETSc to support combining multilevel methods with equation splitting and multirate structure, as well as recent advances in robustness and generality, performance implications of an extremely low-communication variant of multigrid, and a new approach to resilience that enables local reconstruction
of fine-grid state from arbitrarily coarse checkpoints.