Magnetohydrodynamics (MHD) is a fluid theory that describes Plasma Physics by treating the plasma as a fluid of charged particles.  Hence, the equations that describe the plasma form a nonlinear system that couples Navier-Stokes with Maxwell's equations. We describe how the FOSLS method can be applied to incompressible resistive MHD to yield a well-posed, H$^1$-equivalent functional minimization. To solve this system of PDEs, a nested-iteration-Newton-FOSLS-AMG-LAR approach is taken. Much of the work is done on relatively coarse grids, including most of the linearizations.  We show that at most one Newton step and a few V-cycles are all that is needed on the finest grid. Estimates of the local error and of relevant problem parameters that are established while ascending through the sequence of nested grids are used to direct local adaptive mesh refinement (LAR), with the goal of obtaining an optimal grid at a minimal computational cost. An algebraic multigrid solver is used to solve the linearization steps. A parallel implementation is described that uses a binning strategy. We show that once the solution is sufficiently resolved, refinement becomes uniform which essentially eliminates load balancing on the finest grids. The ultimate goal is to resolve as much physics as possible with the least amount of computational work. We show that this is achieved in the equivalent of a few dozen work units on the finest grid. (A work unit equals a fine grid residual evaluation). Numerical results are presented for two instabilities in a large aspect-ratio tokamak, the tearing mode and the island coalescence mode.